aainstitutetext: Department of Physics, Syracuse University, Syracuse, NY
13244, USAbbinstitutetext: Department of Physics and Astronomy, University of
Oklahoma, Norman, OK, 73019, USAccinstitutetext: Department of Physics,
University of Wisconsin, Madison, WI, 53706, USA

# Multi-component Dark Matter and Small Scale Structure Formation

Robert Wiley Deal a Kishan Sankharva b Kuver Sinha a Scott Watson
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider the evolution of non-thermal dark matter perturbations in models
which contain both Weakly Interacting Massive Particles (WIMPs) and axions.
Using constraints from existing observations we examine the percentage of
WIMPs and axions that may comprise the cosmological dark matter budget in
models with an Early Matter Dominated Epoch (EMDE) – where entropy production
is important. After carefully tracking the thermal evolution of the various
species by solving the Boltzmann equations, we consider the enhancement of
perturbations that may have led to early structure formation for axions and
WIMPs. We investigate the impact of enhanced perturbations on the parameter
space of both species, after imposing existing constraints from indirect
detection experiments. Given these constraints we establish the feasibility of
axions to form miniclusters in the early universe in EMDEs for a given
percentage of allowed WIMPs. We find that EMDEs with low reheat temperatures
near the BBN limit are preferred for axion minicluster formation. When the
EMDE is caused by string moduli, the WIMP contribution to the relic density is
set by the moduli branching to dark matter at the level of
$\lesssim\mathcal{O}(10^{-4})$.

##  1 Introduction

The Electroweak hierarchy and Strong CP problems both motivate the presence of
physics Beyond the Standard Model (BSM). In the former case, constructing
particle-based models that address the stabilization of the hierarchy
naturally leads to the prediction of Weakly-Interacting Massive Particles
(WIMPs) giving a candidate for all or part of cosmological Dark Matter (DM).
Whereas, the latter problem also leads to a possible DM candidate – the QCD
axion Weinberg:1977ma ; Peccei:1977hh ; Wilczek:1977pj ; Kim:1979if ;
Shifman:1979if ; Zhitnitsky:1980tq ; Dine:1982ah .

The even more difficult challenge of addressing the origin of the Electroweak
hierarchy requires accounting for the presence of gravity in the theory. This,
along with accounting for adequate cosmic inflation in the early universe
typically implies additional scalar degrees of freedom and additional axion-
like particles (ALPs). Some of these fields can become populated following
inflation and dominate the energy density resulting in a departure from a
strictly thermal history prior to Big Bang Nucleosynthesis (BBN).

Originally posed as the cosmological moduli problem Hall:1983iz ;
deCarlos:1993wie ; Coughlan:1983ci ; Banks:1993en , in recent years these
Early Matter Dominated Epochs (EMDEs) have been viewed more as an opportunity
to probe the early history of the Universe and make connections to string
phenomenology. The primordial origin, abundance, and clustering of WIMPs and
axions can be altered by the pre-BBN cosmic evolution. Moreover, when these
particles decay they can lead to significant entropy production altering the
primordial abundance of WIMPs and axions. EMDEs have been extensively studied
in recent years in the context of their effect on DM searches (for recent
reviews, we refer to Allahverdi:2020bys ; Kane:2015jia ; for earlier studies,
we refer to Giudice:2000ex ; Gelmini:2006pw ) and their connections to the
physics of moduli stabilization Allahverdi:2013noa ; Allahverdi:2014ppa ;
Cicoli:2022uqa ; Dutta:2009uf ; Acharya:2008bk ; Acharya:2009zt .

Much of the previous work on EMDEs has focused on the effects of such eras on
dark sectors consisting of a single component, often a standard WIMP. For
example, the effects on DM indirect and direct detection Blanco:2019eij ;
Easther:2013nga , connections with baryogenesis Kane:2011ih ;
Allahverdi:2010rh ; Allahverdi:2010im ; Allahverdi:2012gk  and structure
formation Miller:2019pss ; Delos:2018ueo ; Delos:2017thv ; Fan:2014zua  have
been probed in detail, yielding a rich phenomenological probe of the early
Universe.

The purpose of this paper is to initiate a study of EMDEs for multicomponent
DM, with an emphasis on the possibility of DM substructure formation. The
focus on enhanced DM sub-structure formation comes from the fact that this is
one of the most important and universal effects of EMDEs: subhorizon DM
perturbations grow linearly during EMDEs and can lead to an enhancement in the
growth of sub-structure on small scales. For a dark sector model with multiple
interacting species or a single self-interacting one, a study combining
alternative cosmological histories with structure formation can be expected to
be quite challenging. As a first step, it is therefore reasonable to restrict
ourselves to models where the different components contributing to the relic
density do not interact with each other. Given the fact that EMDEs typically
arise in string compactifications with both scalar moduli as well as their
pseudoscalar axion-like partners, it is natural to consider the multicomponent
sector to be comprised of WIMPs and axions, specifically the QCD axion. The
phenomenology of such scenaros has been studied in detail, especially in a
supersymmetric setting Baer:2020kwz ; Baer:2019xww ; Bae:2013hma ; Baer:2011uz
; Bae:2014rfa .

Given a multi-component DM sector consisting of WIMPs and axions, one can
focus on the question of sub-structure formation in either sector, and ask
whether the parameter space for such sub-structure formation is consistent
with the observed relic density and various observational bounds. For
concreteness, we will focus on the question of sub-structure formation in the
axion sector (axion miniclusters)111Sub-structure formation in the WIMP sector
has been studied in detail by a subset of the current authors Fan:2014zua .
and study how the existence of such miniclusters affects the properties and
thermal history of the WIMP sector. EMDEs are expected to play a critical role
in the formation of axion miniclusters in the case where Peccei-Quinn (PQ)
symmetry breaking occurs before inflation, which is the case we study. We
first carefully track the thermal history of the modulus, axions, WIMPs and
radiation by solving the coupled Boltzmann equations, and then study the
evolution of the perturbation of modes that enter the Hubble radius during the
EMDE and subsequently through reheating. Imposing the condition for the
formation of miniclusters leads to constraints on the decay temperature of the
modulus, which in turn impact the parameter space available for the WIMP
sector. In particular, we study the resulting constraints on WIMPs coming from
indirect detection data.

The question of axion minicluster formation with EMDEs has garnered some
attention in recent years Nelson:2018via ; Blinov:2019jqc ; Visinelli:2018wza
. Our work differs from these studies mainly by the introduction of the second
component of DM (WIMPs) and by a careful treatment of the interplay between
the components and their thermal histories. As we have argued, string theory
naturally comes with situations where moduli and axions co-exist, and WIMPs
are motivated by a solution to the hierarchy problem. But one can also ask: if
one is interested solely in axion sub-structure with EMDEs, is one forced to
introduce a second DM component? The answer we give is to the affirmative, at
least for a large range of values of the axion decay constant. A careful
computation of the interplay between the axion contribution to the relic
density and the possibility of minicluster formation reveals that sub-
structure formation prefers low modulus reheat temperatures
$T^{\phi}_{D}\lesssim\mathcal{O}(50)$ MeV, which can only accommodate pure
axion DM for decay constants $f_{a}\gtrsim 10^{15}$ GeV. Smaller values of the
decay constant will require a second DM component. These results are depicted
in Fig. 5, and we find that while tuning the initial misalignment can relax
these requirements somewhat, the qualitative trends remain unaltered (a mildly
tuned case is depicted in Fig. 6). The preference of low modulus decay
temperatures for minicluster formation is due to the fact that a prolonged
EMDE is required for perturbations to grow to the requisite level. And it is
precisely in the low reheat regime that axion decay constants below the GUT
scale of $\sim 10^{15}$ GeV will necessitate a second DM component.

In the low reheat regime to which we are forced, the WIMPs have to be produced
from the modulus with small branching ratios $\lesssim\mathcal{O}(10^{-4})$,
corresponding to the “branching scenario” Allahverdi:2020bys ; Kane:2015jia .
For larger branching ratios, WIMP re-annihilation after production from
modulus decay becomes the predominant mechanism for setting the relic density
(the “annihilation scenario”) and such cases are highly constrained by
indirect detection data. Whether the small branching ratios of moduli into
WIMPs required for the branching scenario to be successful can be obtained
from actual string compactificatons is a question we leave for the future.

The paper is structured as follows. In Section 2 we introduce EMDEs and
axions, set our notation, and solve the Boltzmann equations for the various
components. Section 3 addresses substructure formation in the presence of both
axions and WIMPs. We then address the formation of axions in the formation of
miniclusters and then give our results for WIMP constraints in Secction 4. We
then conclude.

##  2 Non-thermal dark matter production

In this section, we discuss the production of both axion and WIMP dark matter
in the context of EMDEs. We begin with a brief review of the QCD axion and its
background evolution, as well as its connections to inflation. We then present
the zeroth-order Boltzmann equations which govern the evolution of the
densities of WIMPs and axions, in addition to radiation and the modulus. This
allows us to analyze the qualitative behavior of the thermal evolution of the
different species and estimate the relic densities in limiting cases. Finally,
we detail our procedure for obtaining numerical solutions of the Boltzmann
equations which we utilize in Sec. 4 and Sec. 5.

###  2.1 The QCD Axion

The QCD axion is a pseudo Nambu-Goldstone boson generated by spontaneous
breaking of the Peccei-Quinn (PQ) symmetry which involves a complex scalar
field called the PQ field. The Lagrangian for the PQ field $\xi$ is given by

$\mathcal{L}_{\xi}=\frac{1}{2}\partial_{\mu}\xi^{\dagger}\,\partial^{\mu}\xi-\frac{\lambda}{4}(\xi^{\dagger}\xi-
f_{a}^{2})^{2},$ (2.1)

where $f_{a}$ is the vacuum expectation value of the PQ field. The expectation
value of the $\xi$ can be written as
$\langle\xi(x)\rangle=f_{a}e^{i\varphi(x)/f_{a}}$ where the phase field
$\varphi(x)$ is identified as the axion field. Here, the axion field
$\varphi(x)$ does not have any potential due to PQ symmetry and therefore, it
can take any value. However, after PQ symmetry breaking, the axion receives a
potential via its coupling to gluons. This potential is minimized at a strong
CP conserving value. Since before PQ transition the axion field takes on any
random value, after the transition its initial value may be ‘misaligned’ from
the minima of its potential. The energy stored in this misaligned axion later
takes the form of dark matter.

The axion mass is expected to be temperature-dependent due to finite-
temperature contributions from instantons, and takes the schematic form222
Here, we have set the domain wall number $N_{\text{DW}}=1$ for simplicity.
This will not have significant impact on our results since most axion models
have $N_{\text{DW}}\sim\mathcal{O}(1-10)$.  Hertzberg:2008wr ; Gross:1980br

$m_{a}(T)=\left(6.2\times 10^{-3}\text{ GeV}\right)\left(\frac{1\text{
GeV}}{f_{a}}\right)\times\begin{cases}1&(T\lesssim\Lambda_{\text{QCD}})\\\
b\left(\frac{\Lambda_{\text{QCD}}}{T}\right)^{4}&(T\gtrsim\Lambda_{\text{QCD}})\end{cases}$
(2.2)

where $b=0.018$ and $\Lambda_{\text{QCD}}\sim 200$ MeV. The axion potential
after PQ symmetry breaking can then be written approximately as

$V_{a}(x,T)=f_{a}^{2}m_{a}^{2}(T)\bigg{[}1-\cos\bigg{(}\frac{\varphi(x)}{f_{a}}\bigg{)}\bigg{]}\simeq\frac{1}{2}m_{a}^{2}(T)\varphi^{2}(x)$
(2.3)

where $m_{a}(T)$ is given by Eq. (2.2).

The background equation of motion for the axion field in an expanding universe
is thus given by

$\ddot{\varphi}+3H(t)\,\dot{\varphi}+m_{a}^{2}(T)\varphi=0$ (2.4)

where $H(t)$ is the Hubble scale, and we use overdots to represent derivatives
with respect to comoving time $t$. Since the axion is minimally coupled, the
above equation can be used to obtain the evolution of the axion energy
density, irrespective of other components in the universe. However, solving
this equation is non-trivial for several reasons. First, note that Eq. (2.4)
is an equation of motion for a damped harmonic oscillator. However, both its
damping coefficient and its mass are effectively time-dependent. This makes a
general analytical solution implausible. Second, the axion oscillates very
rapidly once the damping coefficient ($\sim H$) is smaller than the axion mass
($m_{a}$). This makes Eq. (2.4) a stiff differential equation which is
difficult to solve numerically. Fortunately, we can solve this equation
approximately and deduce the general trend of the evolution of the axion
energy density.

Note that the second term, $3H\dot{\varphi}$, behaves like a time-dependent
friction term which does not allow the axion field to oscillate until $H\simeq
m_{a}$. This gives the condition for critical time $t_{\rm osc}$ at which the
axion oscillations begin:

$H(t_{\rm osc})\simeq m_{a}(T(t_{\rm osc})).$ (2.5)

We can thus take the axion field to be frozen for $t<t_{\rm osc}$, which
produces an equation of state $w_{a}=-1$. The axion field then oscillates very
rapidly for $t>t_{\rm osc}$ such that its average equation of state in this
regime is $\langle w_{a}\rangle=0$. In other words, the axion field behaves
similarly to a cosmological constant before $t_{\rm osc}$ with its number
density $n_{a}=$ constant333 Due to the temperature dependence of the axion
mass, the axion energy density while it is frozen should actually increase as
the universe expands due to its increasing mass. , and it behaves like dark
matter after $t_{\rm osc}$ where its number density dilutes as $n_{a}\sim
a^{-3}$ where we use $a$ to denote the scale factor. Thus, we can estimate the
number density of axions as

$n_{a}(t)\simeq\frac{1}{2}m_{a}(t_{\rm
osc})\,f_{a}^{2}\theta_{i}^{2}\left(\frac{a(t_{\rm osc})}{a(t)}\right)^{3},$
(2.6)

where $\theta_{i}=\varphi_{0}/f_{a}\in[0,\pi)$ is the initial misalignment
angle, which is random and is typically assumed to be $\mathcal{O}(1)$. Now,
the axion energy density at time $t$ is simply given by
$\rho_{a}(t)=m_{a}(t)n_{a}(t)$.

It is worth mentioning that much of this discussion can also apply to the
modulus as well. During the inflationary epoch, moduli are expected to be
displaced from their minima Dine:1995kz ; Dine:1995uk  and frozen in place by
Hubble friction. Once the universe has expanded sufficiently (i.e. when $H\sim
m_{\phi}$), they begin to oscillate and behave as cold matter. However, moduli
typically are not protected by such symmetries and thus are expected to have a
Planck-scale initial amplitude Dine:1995kz ; Dine:1995uk ; Cicoli:2016olq .
Additionally, moduli are typically assumed to be very massive with
$m_{\phi}\gtrsim\mathcal{O}(10-50\text{ TeV})$ to avoid conflicts with BBN.

###  2.2 PQ Breaking Scale and Inflation

The relationship between the PQ breaking scale and the scale of inflation has
been studied both in the context of standard as well as modified cosmologies
Hertzberg:2008wr ; Poulin:2018dzj ; Visinelli:2009kt ; Baer:2011uz ;
Visinelli:2018wza ; Hlozek:2017zzf . Denoting the Hubble scale during
inflation by $H_{I}$, the two scenarios are $f_{a}<H_{I}/2\pi$ (axion absent
during inflation) and $f_{a}>H_{I}/2\pi$ (axion present during inflation). For
the former case, PQ symmetry breaking occurs after inflation and the initial
misalignment angle $\theta_{i}$ acquires different values within the same
Hubble patch. The formation of axion miniclusters in this case has been
studied extensively Kolb:1993zz ; Hogan:1988mp ; Turner:1983sj  and will not
be probed further in the current work.

For the second case, PQ symmetry breaking occurs before inflation, and the
initial misalignment angle takes only a single value in the Hubble horizon. In
this case, axion isocurvature fluctuations are constrained by the cosmic
microwave background and BBN. For standard cosmological histories, the
constraint can be expressed as a relation between $H_{I},\theta_{i},$ and
$f_{a}$ as follows Visinelli:2009zm ; Planck:2018jri :

$H_{I}\,\lesssim\,10^{-5}\,\left(\frac{\Omega_{\rm
CDM}}{\Omega_{a}}\right)\theta_{i}f_{a}\,\,.$ (2.7)

Constraints coming from isocurvature fluctuations in the case of scenarios
with decaying scalar fields were worked out in Visinelli:2009kt . In this
case, one has

$H_{I}\,\lesssim\,1.0\times 10^{8}\left(\frac{\Omega_{\rm
CDM}}{\Omega_{a}}\right)\left(\frac{f_{a}}{10^{12}\,{\rm GeV}}\right)\,\,{\rm
GeV}\,\,,$ (2.8)

which holds for $f_{a}<10^{15}\left(\frac{T^{\phi}_{D}}{1{\rm
MeV}}\right)^{-4/3}$ GeV. For higher values of $f_{a}\gtrsim
10^{15}\left(\frac{T^{\phi}_{D}}{1{\rm MeV}}\right)^{2}$ GeV, one has

$H_{I}\,\lesssim\,5.0\times 10^{10}\left(\frac{\Omega_{\rm
CDM}}{\Omega_{a}}\right)\left(\frac{T^{\phi}_{D}}{1{\rm
MeV}}\right)^{-1/2}\,\,{\rm GeV}\,\,.$ (2.9)

For details of other regimes of $f_{a}$, we refer to Visinelli:2009kt .
Independent of axion physics, the scale of inflation is constrained by Planck
Planck:2018jri  to be

$H_{I}\,\lesssim\,2.5\times 10^{-5}\,m_{P}\,\,.$ (2.10)

For the remainder of this work, we will be agnostic about the details of the
inflationary sector, and assume a value of $H_{I}$ such that all constraints
are satisfied.

###  2.3 The Boltzmann equations

We now shift to the discussion of the zeroth-order Boltzmann equations to
study the production of WIMP and axion dark matter in modulus-dominated non-
thermal cosmologies. We begin with the Friedman, Lemaitre, Robertson, Walker
(FLRW) metric

$ds^{2}=-dt^{2}+a^{2}(t)d\vec{x}^{2},$ (2.11)

with the resulting Hubble equation

$\displaystyle 3H^{2}m_{p}^{2}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\rho_{(\alpha)},$ (2.12)

where $\alpha\in\\{\phi,\chi,r,a\\}$ represents moduli, WIMPs, radiation, and
axions, respectively. We work in units of the reduced Planck mass
$m_{p}=2.44\times 10^{18}~{}\mbox{GeV}$.

For this work, we adopt a coupled Boltzmann equation approach which includes
tracking the particles with all other degrees of freedom treated as radiation.
Instead of working with comoving time, it is convenient to express the full
set of equations as a function of e-foldings, $Hdt=dN=d(\ln{a})$, producing
the set of equations

$\displaystyle\frac{d\rho_{\phi}}{dN}$ $\displaystyle=$
$\displaystyle-3(1+w_{\phi})\rho_{\phi}-\frac{\Gamma_{\phi}}{H}\rho_{\phi},$
(2.13) $\displaystyle\frac{d\rho_{\chi}}{dN}$ $\displaystyle=$
$\displaystyle-3(1+w_{\chi})\rho_{\chi}-\frac{\langle\sigma
v\rangle_{\chi}}{m_{\chi}H}\left(\rho_{\chi}^{2}-\overline{\rho}_{\chi}^{2}\right)+B_{\phi\rightarrow\chi}\left(\frac{\Gamma_{\phi}}{H}\right)\rho_{\phi},$
(2.14) $\displaystyle\frac{dn_{a}}{dN}$ $\displaystyle=$
$\displaystyle-3\left(1+w_{a}\right)n_{a},$ (2.15)
$\displaystyle\frac{d\rho_{r}}{dN}$ $\displaystyle=$
$\displaystyle-4\rho_{r}+\frac{\langle\sigma
v\rangle_{\chi}}{m_{\chi}H}\left(\rho_{\chi}^{2}-\overline{\rho}_{\chi}^{2}\right)+(1-B_{\phi\rightarrow\chi})\frac{\Gamma_{\phi}}{H}\rho_{\phi},$
(2.16)

subject to the energy constraint Eq. (2.12) and where we have kept the
equation of state $w_{\alpha}$ general. Note that we have also expressed the
equation governing the evolution of the axion in terms of its number density,
$n_{a}$, which is simpler to work with than its energy density.444 The
Boltzmann equation governing the evolution of the axion’s energy density
requires the inclusion of a term of the form $dm_{a}/dN$ due to its
aforementioned temperature-dependent mass.  Indeed, because the axion field is
always non-relativistic we can always retrieve its energy density from
knowledge of its number density and its mass. Before we continue, it is worth
discussing the limiting cases of the above Boltzmann equations which reduce to
familiar expressions typically used in the literature.

###  2.4 Qualitative behavior and limiting cases

Our first task is to estimate the scale of modulus decay. From the Boltzmann
equation of the modulus, Eq. (2.13), we see the decay term becomes sizeable
when

$\Gamma_{\phi}\sim H.$ (2.17)

In the sudden decay approximation and assuming a radiation-dominated universe
after the modulus decays, we can estimate the decay temperature $T_{D}^{\phi}$
as

$T_{D}^{\phi}\simeq\sqrt{m_{P}\Gamma_{\phi}}\left(\frac{90}{\pi^{2}g_{*}(T_{D}^{\phi})}\right)^{1/4}.$
(2.18)

This temperature is also frequently referred to as the “reheating temperature”
in much of the literature. As moduli are expected to be only gravitationally-
coupled, we note that the decay temperature can be expected to be quite low.
If $T_{D}^{\phi}\lesssim\mathcal{O}(1\text{ MeV})$, decay occurs during or
after Big Bang Nucleosynthesis (BBN) Kawasaki:1999na ; Kawasaki:2000en ;
Hasegawa:2019jsa  and - due to the substantial entropy production - erases the
successful predictions of BBN. This is the standard cosmological moduli
problem (CMP) and sets a lower bound on the decay width $\Gamma_{\phi}\gtrsim
5\times 10^{-25}$ GeV, which in turn sets a constraint on the modulus mass for
a particular set of couplings. As a precise limit on the modulus decay
temperature is model-dependent Hasegawa:2019jsa , we adopt a naive BBN cutoff
of 1 MeV throughout.

It is also useful to estimate the scale at which entropy begins to be injected
into the thermal bath. To proceed, we first define the temperature at which
the modulus begins to dominate the energy density of the universe as
$T_{e}^{\phi}$. From imposing $\rho_{r}\sim\rho_{\phi}$ and assuming a
radiation-dominated universe before modulus domination, the temperature of
modulus-radiation equality can be shown to be555 This form assumes that
modulus oscillations onset only after any inflationary reheating has finished
so that the modulus is unaffected by entropy production from the decay of the
inflaton. See e.g. Bae:2022okh ; WileyDeal:2023sry  for details.  Bae:2022okh

$T_{e}^{\phi}\simeq\frac{3}{2}\left(\frac{10}{\pi^{2}g_{*}(T_{\text{osc}}^{\phi})}\right)^{1/4}\left(\frac{\phi_{0}}{m_{P}}\right)^{2}\sqrt{m_{\phi}m_{P}}$
(2.19)

where $\phi_{0}$ is the initial amplitude of the modulus. As the universe is
matter-dominated between $T_{e}^{\phi}$ and $T_{D}^{\phi}$, it is a simple
matter to relate the Hubble scales at both eras. Upon noting that entropy is
conserved between $T_{e}^{\phi}$ and - by definition - the entropy injection
temperature $T_{S}^{\phi}$, while entropy is not conserved between
$T_{S}^{\phi}$ and $T_{D}^{\phi}$, we can arrive at the following estimate for
the entropy injection temperature:

$T_{S}^{\phi}\simeq\left(T_{e}^{\phi}\left(T_{D}^{\phi}\right)^{4}\right)^{1/5}.$
(2.20)

Although this entropy injection temperature will not be necessary for our
discussion on WIMP production, it will become a key but implicit feature in
determining the axion relic density within a modulus-dominated cosmology.

We now consider the limiting cases for WIMP production, beginning with the
case where WIMP annihilations are dominant. Considering the regime close to
(or below) the thermal freeze-out scale and assuming the WIMPs behave as cold
dark matter in this regime (i.e. $w_{\chi}=0$ and $\rho_{\chi}\simeq
m_{\chi}n_{\chi}$), the Boltzmann equation for WIMPs can be rewritten in terms
of the abundance yield $Y_{\chi}\equiv n_{\chi}/s$ as

$\frac{dY_{\chi}}{dN}\simeq-\frac{\langle\sigma
v\rangle_{\chi}}{H}Y_{\chi}^{2}s.$ (2.21)

Taking the limit of constant $\langle\sigma v\rangle_{\chi}$, this equation
can be straightforwardly integrated from the scale of modulus decay
$T_{D}^{\phi}$ to the present time to give the result

$Y_{\chi}\simeq\frac{H(T_{D}^{\phi})}{\langle\sigma
v\rangle_{\chi}\,s(T_{D}^{\phi})}.$ (2.22)

Assuming no significant source of entropy production after the modulus has
decayed, the WIMP relic density can be estimated by Kolb:1990vq

$\Omega_{\chi}h^{2}\simeq\frac{m_{\chi}Y_{\chi}h^{2}}{\rho_{c}/s_{0}}$ (2.23)

where $s_{0}\sim 2970\text{ cm}^{-3}$ and $\rho_{c}/h^{2}\sim 1.05\times
10^{-5}\text{ GeV}\cdot\text{cm}^{-3}$. This can also be used to relate the
produced WIMP relic density $\Omega_{\chi}h^{2}$ to the expected WIMP relic
density produced in a standard thermal scenario
$\Omega_{\chi}^{\text{th}}h^{2}$:

$\Omega_{\chi}h^{2}=\Omega_{\chi}^{\text{th}}h^{2}\times\max\\{T_{f}/T_{D}^{\phi},1\\}.$
(2.24)

This well-known result is typically coined the annihilation scenario in the
literature Allahverdi:2013noa ; Allahverdi:2014ppa ; Cicoli:2022uqa .

The other limiting case for the production of WIMPs is when the modulus decay
term dominates. Making again the assumption that the WIMPs behave as cold dark
matter, the Boltzmann equation for WIMPs can be written as

$\frac{dY_{\chi}}{dN}=B_{\phi\rightarrow\chi}\frac{\Gamma_{\phi}}{H}Y_{\phi}$
(2.25)

where we again rewrite the Boltzmann equation in terms of the abundance yield
$Y_{\chi}$. Integrating this equation from the scale of modulus decay to the
present time then gives the result

$Y_{\chi}\simeq B_{\phi\rightarrow\chi}Y_{\phi}(T_{D}^{\phi}).$ (2.26)

Once again, the WIMP relic density can be estimated with Eq. (2.23). This
result is typically denoted as the branching scenario Allahverdi:2013noa ;
Allahverdi:2014ppa ; Cicoli:2022uqa .

We now turn to a discussion of the axion. In this work, we concern ourselves
only with the case where the modulus decays after the axion field has begun to
oscillate, i.e. $T_{D}^{\phi}<T_{\text{osc}}^{a}$, so that axion miniclusters
might form if PQ symmetry breaking occurs at or above the inflationary scale
Nelson:2018via . This has the crucial property that entropy dilution from the
modulus will reduce the late-time relic density of the misalignment axions
Baer:2023bbn . In a “decay-dominated” background - i.e. during the era in
which the modulus injects significant entropy into the thermal bath - we can
then estimate the Hubble scale at the onset of axion oscillations in terms of
the modulus decay temperature as Kolb:1990vq

$H(T_{\text{osc}}^{a})\simeq\frac{g_{*}(T_{\text{osc}}^{a})}{g_{*}(T_{D}^{\phi})}\frac{(T_{\text{osc}}^{a})^{4}}{(T_{D}^{\phi})^{2}}\sqrt{\frac{\pi^{2}g_{*}(T_{D}^{\phi})}{90m_{P}^{2}}}$
(2.27)

where we assume a transition to a radiation-dominated universe at
$T_{D}^{\phi}$. It is then straightforward to compute $T_{\text{osc}}^{a}$
from Eq. (2.2) and Eq. (2.5), which we find to be roughly
$T_{\text{osc}}^{a}\sim\mathcal{O}(0.1-1\text{ GeV})$ depending on the value
of $T_{D}^{\phi}$. Having found the axion oscillation temperature, we can then
estimate the axion number density when oscillations begin by

$n_{a}(T_{\text{osc}}^{a})\simeq\frac{1}{2}m_{a}(T_{\text{osc}}^{a})\mathcal{A}_{0}^{2}$
(2.28)

where, accounting for anharmonic effects, the initial oscillation amplitude
$\mathcal{A}_{0}$ is given by Visinelli:2009zm

$\mathcal{A}_{0}^{2}=1.44\,f_{a}^{2}\theta_{i}^{2}\,\left[\log\left(\frac{e}{1-\theta_{i}^{2}/\pi^{2}}\right)\right]^{7/6}$
(2.29)

in terms of the Peccei-Quinn scale $f_{a}$ and the initial misalignment angle
$\theta_{i}$. Accounting for the remaining matter-dominated period between
$T_{\text{osc}}^{a}$ and $T_{D}^{\phi}$ so that

$n_{a}(T_{D}^{\phi})=n_{a}(T_{\text{osc}}^{a})\left(\frac{a(T_{\text{osc}}^{a})}{a(T_{D}^{\phi})}\right)^{3}\simeq
n_{a}(T_{\text{osc}}^{a})\left(\frac{g_{*}(T_{D}^{\phi})(T_{D}^{\phi})^{4}}{g_{*}(T_{\text{osc}}^{a})(T_{\text{osc}}^{a})^{4}}\right)^{2}$
(2.30)

and assuming no entropy production after $T_{D}^{\phi}$, we arrive at the
relic density estimate for misalignment axions:

$\Omega_{a}h^{2}\simeq\frac{45}{4\pi^{2}}\frac{m_{a}(T=0)\,m_{a}(T_{\text{osc}}^{a})\,\mathcal{A}_{0}^{2}}{\rho_{c}/(h^{2}s_{0})}\frac{g_{*}(T_{D}^{\phi})}{g_{*}^{2}(T_{\text{osc}}^{a})}\frac{(T_{D}^{\phi})^{5}}{(T_{\text{osc}}^{a})^{8}}.$
(2.31)

###  2.5 Numerical solutions to Boltzmann equations

Now that we have a few qualitative estimates for the scales and predicted
abundances, we can return to the procedure of numerically integrating Eqs.
(2.13)-(2.16). A small codebase was written in Python incorporating these
differential equations, which were solved using a backward differentiation
formula (BDF) implementation in SciPy’s ode class 2020SciPy-NMeth . At each
step, we recover the temperature from the radiation energy density. It is then
a simple matter to recover the current step’s Hubble parameter and the axion
mass, as well as the equilibrium energy density of the WIMP which, making the
conventional definition $x\equiv m_{\chi}/T_{\chi}$, can be approximated as
Kolb:1990vq ; Baer:2011uz

$\overline{\rho}_{\chi}=\begin{cases}\frac{7}{8}\frac{\pi^{2}}{30}\,gT^{4}&(3/2\leq
x^{-1})\\\
\frac{1}{2\pi^{2}}\,gT^{4}\,x^{2}\left(xK_{1}(x)+3K_{2}(x)\right)&(1/10<x^{-1}<3/2)\\\
\frac{1}{(2\pi)^{3/2}}\,gT^{4}\,x^{5/2}\exp(-x)&(x^{-1}\leq 1/10)\end{cases}$
(2.32)

where the middle case interpolates between the relativistic and non-
relativistic limits, and where $g=2$ for a spin-$1/2$ WIMP.

Within the Boltzmann equation implementation, our code adopts the equation of
state

$w_{\alpha}=\begin{cases}-1&(m_{\alpha}<H)\\\ 0&(m_{\alpha}\geq H)\end{cases}$
(2.33)

for both the modulus and axion (i.e. $\alpha\in\\{\phi,a\\}$), where the axion
uses its mass for the current step. The equation of state for WIMPs, however,
is slightly more complicated as we numerically integrate from the regime where
WIMPs are relativistic to the regime where they become cold. The WIMP equation
of state can be approximated as Kolb:1990vq ; Baer:2011uz

$w_{\chi}=\begin{cases}1/3&(3/2\leq x^{-1})\\\
\left(x\frac{K_{1}(x)}{K_{2}(x)}+3\right)^{-1}&(1/10<x^{-1}<3/2)\\\
0&(x^{-1}\leq 1/10)\end{cases}$ (2.34)

where we again use $x\equiv m_{\chi}/T_{\chi}$. However, there is an
additional complication - when WIMPs are produced via modulus decay, they may
be highly relativistic even at temperatures below the standard thermal freeze-
out temperature simply due to the comparatively large modulus mass
WileyDeal:2023sry . Thus, the WIMP equation of state requires a subtle
distinction in its temperature. In this work, we accomplish this by
incorporating an additional Boltzmann equation governing the number density of
the WIMP. This additional equation is entirely independent of the equation of
state, which allows us to both numerically fit $T_{\chi}$ through a comparison
of $\rho_{\chi}$ to $n_{\chi}m_{\chi}$ as well as provide a secondary estimate
for the WIMP relic density. Interestingly, we find agreement between the
number density and energy density Boltzmann equations to well below the
percent-level for the calculated WIMP relic density.

All that remains is to determine the initial conditions. As the results are
expected to be independent of the inflationary reheating temperature666 This
behavior was confirmed in Baer:2023bbn  for predictions of both WIMP and axion
dark matter, as well as production of dark radiation, using a similar
codebase.  $T_{R}$, we set our initial conditions at this point. Assuming that
the universe is radiation-dominated at $T_{R}$ (with a period of EMD driven by
the modulus occurring shortly after), the initial conditions follow
immediately. The initial radiation energy density is set simply by

$\rho_{r}^{(0)}=\frac{\pi^{2}}{30}g_{*}(T_{R})T_{R}^{4}.$ (2.35)

The initial energy density of the modulus is given by

$\rho_{\phi}^{(0)}=\frac{1}{2}m_{\phi}^{2}\phi_{0}^{2}$ (2.36)

where we take $\phi_{0}=m_{P}$, as expected on dimensional grounds Dine:1995kz
; Dine:1995uk  and was found in an explicit calculation for Kähler inflation
Cicoli:2016olq . We then take Eq. (2.28) for the initial number density of the
axion - which is updated at each step until oscillations begin to account for
the varying mass. Finally, a quick comparison of the interaction rate
$\langle\sigma v\rangle_{\chi}\overline{\rho}_{\chi}(T_{R})/m_{\chi}$ to the
Hubble rate $H(T_{R})$ for the WIMP cross sections and masses of interest
suggest that the WIMPs begin in thermal equilibrium. Thus, we take
$\rho_{\chi}^{(0)}=\overline{\rho}_{\chi}(T_{R})$.

Figure 1:  Evolution of energy densities $\rho$ versus number of $e$-foldings
$N$ in a scenario resembling the annihilation scenario for the WIMP. Here, we
take $m_{\phi}=100$ TeV, $m_{\chi}=200$ GeV, $c=27.5$,
$B_{\phi\rightarrow\chi}=0.01$, and $\langle\sigma v\rangle_{\chi}=1.05\times
10^{-25}\text{ cm}^{3}\cdot\text{s}^{-1}$. We also take $f_{a}=10^{11}$ GeV,
$\theta_{i}=3.113$, and $T_{R}=10^{12}$ GeV.

The remaining parameters in Eqs. (2.14)-(2.16) are then taken as free
parameters, yielding the input set

$\\{c,\,m_{\phi},\,B_{\phi\rightarrow\chi},\,\langle\sigma
v\rangle_{\chi},\,m_{\chi},\,f_{a},\,\theta_{i},\,T_{R}\\}$ (2.37)

where we define $c$ through the expected decay rate formula Bae:2022okh ;
Cicoli:2022uqa ; Baer:2023bbn

$\Gamma_{\phi}\simeq\frac{c}{48\pi}\frac{m_{\phi}^{3}}{m_{P}^{2}}.$ (2.38)

The number of relativistic degrees of freedom, $g_{*}(T)$, is fitted from data
used in Baer:2011uz ; Baer:2023bbn .

Figure 2:  Evolution of energy densities $\rho$ versus number of $e$-foldings
$N$ in a scenario resembling the branching scenario for the WIMP. Here, we
take $m_{\phi}=100$ TeV, $m_{\chi}=200$ GeV, $c=27.5$,
$B_{\phi\rightarrow\chi}=10^{-8}$, and $\langle\sigma
v\rangle_{\chi}=1.05\times 10^{-25}\text{ cm}^{3}\cdot\text{s}^{-1}$. We also
take $f_{a}=10^{11}$ GeV, $\theta_{i}=3.113$, and $T_{R}=10^{12}$ GeV.

We display the evolution of energy densities using the numeric Boltzmann
solutions for parameters that predominantly give the annihilation scenario in
Fig. 1 and the branching scenario in Fig. 2. For both figures, we adopt
parameter values motivated by our demands that 1.) the WIMP annihilation cross
section is consistent with experimental data, 2.) an EMD phase is prolonged
sufficiently to allow axion miniclusters to form based on arguments given by
Nelson et al. Nelson:2018via  and Blinov et al. Blinov:2019jqc , 3.) the
modulus decays before BBN, and 4.) the modulus has an effective coupling
parameter $c$ that is compatible with explicit models.777We take $c=27.5$
throughout this work based on results presented in Cicoli:2022uqa . Although
it is likely that $c$ could be much larger by the inclusion of additional
sectors to which the modulus can decay, our chosen value is near what one
might expect to be a lower bound for the effective coupling as it includes
only Standard Model degrees of freedom. As $\Gamma_{\phi}$ is proportional to
$c$, we will see that raising $c$ significantly leads to a larger
$T_{D}^{\phi}$ and thus makes our results even more constrained. We also
display a comparison to the standard thermal evolution (where the modulus is
absent), which are shown by the dashed lines of the respective constituents’
colors. The vertical dashed lines in both figures identify qualitative shifts
in the cosmology.

##  3 Substructure formation

To understand the formation of substructure, we need to consider perturbations
of the background equations Eqs. (2.13) – (2.16). The metric in longitudinal
gauge takes the form

$ds^{2}=-\left(1+2\Phi\right)dt^{2}+a(t)^{2}\left(1-2\Psi\right)\delta_{ij}dx^{i}dx^{j}.$
(3.1)

In the absence of anisotropic stress for the fluid sources, we have
$\Phi=\Psi$. Working in momentum space and suppressing the wave number, the
time-time component of the perturbed Einstein equation is

$\left(\frac{k^{2}}{3a^{2}H^{2}}+1\right)\Phi+\Phi^{\prime}=-\frac{1}{6H^{2}m_{p}^{2}}\sum_{\alpha}\delta\rho_{(\alpha)},$
(3.2)

where the prime denotes derivatives with respect to number of e-folds $N$, and
$v_{(\alpha)},\delta\rho_{(\alpha)}$, $\delta p_{(\alpha)}$ are scalar
velocity, density and pressure perturbations for each fluid, respectively.
Introducing fractional density perturbations
$\delta_{(\alpha)}\equiv\delta\rho_{(\alpha)}/\rho_{(\alpha)}$ and defining
the velocity perturbation for each fluid as
$\theta_{(\alpha)}=a^{-1}\nabla^{2}v_{(\alpha)}$, the continuity equations in
momentum space are given by

$\displaystyle\delta^{\prime}_{\phi}+(1+w_{\phi})\left(\frac{\theta_{\phi}}{aH}-3\Phi^{\prime}\right)$
$\displaystyle=$ $\displaystyle-\frac{\Gamma_{\phi}}{H}\Phi$ (3.3)
$\displaystyle\delta^{\prime}_{\chi}+(1+w_{\chi})\left(\frac{\theta_{\chi}}{aH}-3\Phi^{\prime}\right)$
$\displaystyle=$ $\displaystyle
B_{\phi\rightarrow\chi}\left(\frac{\Gamma_{\phi}}{H}\right)\left(\frac{\rho_{\phi}}{\rho_{\chi}}\right)\left[\delta_{\phi}-\delta_{\chi}+\Phi\right]$
(3.4) $\displaystyle-\frac{\langle\sigma
v\rangle_{\chi}}{m_{\chi}H}\left(\rho_{\chi}\left(\delta_{\chi}+\Phi\right)-\frac{\overline{\rho}^{2}_{\chi}}{\rho_{\chi}}\left(\Phi+2\overline{\delta}_{\chi}-\delta_{\chi}\right)\right)$
$\displaystyle\delta^{\prime}_{a}+(1+w_{a})\left(\frac{\theta_{\chi}}{aH}-3\Phi^{\prime}\right)$
$\displaystyle=$
$\displaystyle\frac{dm_{a}}{dN}\frac{1}{m_{a}}\Phi-3c^{2}_{s,a}\delta_{a}$
(3.5)
$\displaystyle\delta^{\prime}_{r}+\frac{4}{3}\left(\frac{\theta_{r}}{aH}-3\Phi^{\prime}\right)$
$\displaystyle=$
$\displaystyle(1-B_{\phi\rightarrow\chi})\left(\frac{\Gamma_{\phi}}{H}\right)\left(\frac{\rho_{\phi}}{\rho_{r}}\right)\left[\delta_{\phi}-\delta_{r}+\Phi\right]$
(3.6) $\displaystyle+$ $\displaystyle\frac{\langle\sigma
v\rangle_{\chi}}{m_{\chi}H}\left(\frac{\rho_{\chi}^{2}}{\rho_{r}}\right)\left(\left(2\delta_{\chi}-\delta_{r}+\Phi\right)-\frac{\overline{\rho}_{\chi}^{2}}{\rho^{2}_{\chi}}\left(2\delta_{\chi}-\delta_{r}+\Phi\right)\right).\;\;\;\;\;$

Similarly, the equations for velocity perturbations are

$\displaystyle\left(1+w_{\phi}\right)\left(\theta_{\phi}^{\prime}+4\theta_{\phi}-\frac{k^{2}}{aH}\Phi\right)-3\theta_{\phi}$
$\displaystyle=$ $\displaystyle 0$ (3.7)
$\displaystyle\left(1+w_{\chi}\right)\left(\theta_{\chi}^{\prime}+4\theta_{\chi}-\frac{k^{2}}{aH}\Phi\right)-3\theta_{\chi}$
$\displaystyle=$ $\displaystyle
B_{\phi\rightarrow\chi}\left(\frac{\Gamma_{\phi}}{H}\right)\left(\frac{\rho_{\phi}}{\rho_{\chi}}\right)\left[\theta_{\phi}-\theta_{\chi}\right]$
(3.8)
$\displaystyle\left(1+w_{a}\right)\left(\theta_{a}^{\prime}+4\theta_{a}-\frac{k^{2}}{aH}\Phi\right)-3\theta_{a}$
$\displaystyle=$
$\displaystyle\left(3-\frac{dm_{a}}{dN}\frac{1}{m_{a}}\right)c^{2}_{s,a}\theta_{a}+\frac{k^{2}}{aH}c_{s,a}^{2}\delta_{a}\;\;\;\;\;\;$
(3.9)
$\displaystyle\theta_{r}^{\prime}-\frac{k^{2}}{aH}\left(\frac{\delta_{r}}{4}+\Phi\right)$
$\displaystyle=$
$\displaystyle\left(1-B_{\phi\rightarrow\chi}\right)\left(\frac{\Gamma_{\phi}}{H}\right)\left(\frac{\rho_{\phi}}{\rho_{r}}\right)\left[\frac{3}{4}\theta_{\phi}-\theta_{r}\right]\;\;\;\;\;\;$
(3.10) $\displaystyle+\frac{\langle\sigma
v\rangle_{\chi}}{m_{\chi}H}\left(\frac{\rho_{\chi}^{2}-\overline{\rho}_{\chi}^{2}}{\rho_{r}}\right)\left(\frac{3}{4}\theta_{\chi}-\theta_{r}\right).$

This set of differential equations can be closed by the perturbed Einstein
equation Eq. (3.2).

To calculate the evolution of perturbations, we need to specify initial
conditions for the modes. Following Fan:2014zua , we set the modes after
modulus domination when all modes of interest are super-Hubble, $k/aH\to 0$.
As discussed in Fan:2014zua , isocurvature constraints do not apply to the
model we are considering here and so we are interested in strictly adiabatic
initial conditions for the multi-fluid perturbations

$\frac{\delta\rho^{(0)}_{\alpha}}{\rho^{\prime}_{\alpha}}=\frac{\delta\rho^{(0)}_{\beta}}{\rho^{\prime}_{\beta}}.$
(3.11)

We will initially neglect the axion sector as its perturbations decouple from
the system above – we will address the axion perturbations in the next section
when discussing miniclusters. To establish the initial conditions for the
modulus and WIMP sectors, it is useful to find approximate analytic solutions
to background well within the modulus dominated phase with $t\sim H^{-1}\gg
m_{\phi}^{-1}$. An approximate solution is then

$\displaystyle\rho_{\phi}(N)$ $\displaystyle\simeq$
$\displaystyle\rho_{\phi}^{(0)}~{}e^{-3N}$ (3.12)
$\displaystyle\rho_{\chi}(N)$ $\displaystyle=$
$\displaystyle\rho_{\chi}^{(0)}~{}e^{-3N/2}$ (3.13) $\displaystyle\rho_{r}(N)$
$\displaystyle=$ $\displaystyle\rho_{r}^{(0)}~{}e^{-3N/2}$ (3.14)

where we choose initial values so that
$\rho_{\chi}^{(0)},\rho_{r}^{(0)}<<\rho_{\phi}^{(0)}$ and WIMP DM will be
primarily of non-thermal origin (i.e., we neglect equilibrium terms). The
scaling behavior in Eq. (3.13) and Eq. (3.14) results from the near
cancellation between decays and WIMP annihilation whereas the axion sector
decouples and can evolve as matter or radiation.

Using the background fluid equations with our approximation, and remembering
that during modulus domination $\Gamma_{\phi}/H\ll 1$, it follows that

$\delta_{\phi}^{(0)}=2\delta_{\chi}^{(0)}=2\delta_{r}^{(0)}.$ (3.15)

This relation differs from the standard relations due to the presence of
decays and entropy production. Taking the super-horizon limit $k/aH\to 0$ of
Eq. (3.2) in a modulus dominated universe, we have

$\Phi\simeq-\frac{1}{6H^{2}m_{p}^{2}}\rho_{\phi}\delta_{\phi}$ (3.16)

where we used conservation of $\Phi$ on super-Hubble scales. Since
$\rho_{\phi}\simeq 3H^{2}m_{p}^{2}$ during modulus domination, Eq. (3.16)
implies the following initial condition for long wavelength gravitational
perturbations $\delta_{\phi}^{(0)}=-2\Phi_{0}$ and it follows from Eq. (3.11)
that $\delta_{\chi}^{(0)}=\delta_{r}^{(0)}=-\Phi_{0}$. Finally, because scalar
velocity perturbations quickly decay outside of the Hubble radius their
initial values can be ignored.

###  3.1 Perturbation evolution during modulus domination

In this section, we examine the evolution of the perturbations for modes that
enter the Hubble radius during modulus domination. We note that these modes
will be small compared to the size of the horizon at reheating,
$k^{-1}<k_{rh}^{-1}$, and thus it will be important for determining the growth
of structure at that time.

Modulus perturbations

During modulus domination, we have an effectively matter dominated universe
and can therefore set $\Phi=\Phi_{0}$ for both super and sub-Hubble scales.
Using this, we can rewrite Eq. (3.7) as

$\theta_{\phi}^{\prime}+\theta_{\phi}=\frac{k^{2}}{H_{0}}\Phi_{0}e^{N/2},$
(3.17)

where we used $H=H_{0}~{}e^{-3N/2}$ in a matter dominated Universe. This
equation can be solved to give the behavior for all wavelengths. Concentrating
on the growing mode, we have

$\theta_{\phi}(k,N)=\frac{2}{3}\frac{k^{2}}{H_{0}}\Phi_{0}e^{N/2},$ (3.18)

which confirms that long wavelength vector modes are unimportant. From Eq.
(3.18), we can derive the evolution of the modulus perturbation
$\delta_{\phi}$ during the modulus dominated era subject to the initial
condition $\delta_{\phi}^{(0)}=-2\Phi_{0}$. Noting that until the time of
reheating we have $\Gamma_{\phi}/H\ll 1$ and $\Phi$ is constant, we can
rewrite Eq. (3.3) as

$\delta^{\prime}_{\phi}(k,N)\simeq-\frac{\theta_{\phi}}{H_{0}}e^{N/2}$ (3.19)

and using the result in Eq. (3.18), we integrate to find

$\delta_{\phi}(k,N)\simeq-2\Phi_{0}-\frac{2}{3}\frac{k^{2}}{H_{0}^{2}}\Phi_{0}e^{N}$
(3.20)

which is again valid on all scales.

WIMP and radiation perturbations

In the absence of source terms for the perturbation equations, the
perturbations evolve as expected in a matter dominated universe. However, the
additional terms will be important during the period of modulus domination.
Solutions for the complete system can be found as in Fan:2014zua  where by
using the background approximation Eqs. (3.12) - (3.14), we have time-
dependent coefficients scaling as

$\displaystyle
B_{\chi}\frac{\Gamma_{\phi}}{H}\left(\frac{\rho_{\phi}}{\rho_{\chi}}\right)$
$\displaystyle\longrightarrow$ $\displaystyle
B_{\chi}\frac{\Gamma_{\phi}}{H_{0}}\left(\frac{\rho_{\phi}^{(0)}}{\rho_{\chi}^{(0)}}\right)\equiv
A_{1}$ (3.21) $\displaystyle\frac{\langle\sigma
v\rangle}{m_{\chi}H}\rho_{\chi}$ $\displaystyle\longrightarrow$
$\displaystyle\frac{\langle\sigma
v\rangle}{m_{\chi}H_{0}}\rho_{\chi}^{(0)}\equiv A_{2}.$ (3.22)

for density perturbations of the WIMP. For density perturbations of radiation,
from Eq. (3.6) and Eq. (3.10) we again find that the scaling cancels and the
coefficients are determined by their initial values,

$\displaystyle(1-B_{\chi})\frac{\Gamma_{\phi}}{H}\left(\frac{\rho_{\phi}}{\rho_{r}}\right)$
$\displaystyle\longrightarrow$
$\displaystyle(1-B_{\chi})\frac{\Gamma_{\phi}}{H_{0}}\left(\frac{\rho_{\phi}^{(0)}}{\rho_{r}^{(0)}}\right)\equiv
A_{3}$ (3.23) $\displaystyle\frac{\langle\sigma
v\rangle}{m_{\chi}H}\left(\frac{\rho_{\chi}}{\rho_{r}}\right){\rho_{\chi}}$
$\displaystyle\longrightarrow$ $\displaystyle\frac{\langle\sigma
v\rangle}{m_{\chi}H_{0}}\left(\frac{\rho^{(0)}_{\chi}}{\rho_{r}^{(0)}}\right){\rho^{(0)}_{\chi}}\equiv
A_{4}.$ (3.24)

It was shown in Fan:2014zua  that when selecting a range of initial values
motivated from SUSY model building, the annihilation and decay terms are of
comparable importance.

The velocity perturbations of the WIMP fluid during modulus domination can be
found by using Eq. (3.18) in Eq. (3.8), producing the result Fan:2014zua

$\theta_{\chi}(k,N)=\frac{2}{3}\frac{k^{2}}{H_{0}}\Phi_{0}e^{N/2}.$ (3.25)

Similarly, using Eq. (3.20) and Eq. (3.25) in Eq. (3.5) and remembering that
the background coefficients are constants, one finds

$\delta_{\chi}(k,N)=-\Phi_{0}-\frac{2}{3}\left(\frac{1+A_{1}}{1+A_{1}+A_{2}}\right)\frac{k^{2}}{H_{0}^{2}}\Phi_{0}e^{N}$
(3.26)

which again is valid on both super-Hubble and sub-Hubble scales, and we have
used the initial conditions
$\delta_{\phi}^{(0)}=2\delta_{\chi}^{(0)}=-2\Phi_{0}$.

Given the solutions for the modulus and WIMP perturbations, we can solve for
the radiation fluid perturbations. The equation for radiation perturbations is
Fan:2014zua

$\delta_{r}^{\prime\prime}+\left(2A-\frac{1}{2}\right)\delta_{r}^{\prime}+\left(A^{2}-\frac{A}{2}+\frac{k^{2}}{3H_{0}^{2}}e^{N}\right)\delta_{r}=S(N),$
(3.27)

where $A\equiv A_{3}+A_{4}$ and

$\alpha=\frac{2}{3}\left(\frac{A_{1}A_{3}+2A_{1}A_{4}+A_{2}A_{3}+2A_{4}+A_{3}}{1+A_{1}+A_{2}}\right),$
(3.28)

and the source term is given by

$S(N)\equiv-\left(A^{2}-\frac{A}{2}\right)\Phi_{0}-\left(\frac{\alpha}{2}\left(2A+1\right)+\frac{2}{3}\left(A+2\right)\right)\frac{k^{2}}{H_{0}^{2}}\Phi_{0}e^{N}.$
(3.29)

A solution to Eq. (3.27) is

$\frac{\delta_{r}}{\Phi_{0}}=-e^{-AN}\bigg{[}\cos\bigg{\\{}\frac{2}{\sqrt{3}}\frac{k}{H_{0}}(e^{N/2}-1)\bigg{\\}}+\sqrt{3}\frac{AH_{0}}{k}\sin\bigg{\\{}\frac{2}{\sqrt{3}}\frac{k}{H_{0}}(e^{N/2}-1)\bigg{\\}}\bigg{]}\\\
+\sqrt{3}\frac{H_{0}}{k}e^{-AN}(2A^{2}-A)\int_{1}^{e^{N/2}}x^{2A-2}\sin\bigg{\\{}\frac{2}{\sqrt{3}}\frac{k}{H_{0}}(x-e^{N/2})\bigg{\\}}{\rm
d}x\\\
+\frac{1}{\sqrt{3}}\frac{k}{H_{0}}e^{-AN}(6A\alpha+4A+3\alpha+8)\int_{1}^{e^{N/2}}x^{2A}\sin\bigg{\\{}\frac{2}{\sqrt{3}}\frac{k}{H_{0}}(x-e^{N/2})\bigg{\\}}{\rm
d}x,$ (3.30)

which can be further simplified by integration by parts. We thus obtain

$\frac{\delta_{r}}{\Phi_{0}}=-1+\frac{1}{\sqrt{3}}\frac{k}{H_{0}}[(3\alpha+2)(2A+1)+4]\int_{1}^{e^{N/2}}x^{2A}\sin\bigg{\\{}\frac{2}{\sqrt{3}}\frac{k}{H_{0}}(x-e^{N/2})\bigg{\\}}{\rm
d}x.$ (3.31)

In the absence of decays and annihilation (corresponding to $A=\alpha=0$
above), the exact solution to Eq. (3.27) simplifies significantly for all
scales:

$\delta_{r}=-4\Phi_{0}+3\Phi_{0}\cos\left(\frac{2k}{\sqrt{3}H_{0}}(e^{N/2}-1)\right).$
(3.32)

After passing through the Hubble radius, the radiation modes begin to
oscillate with a maximum amplitude reached when the source term in Eq. (3.29)
is in resonance with the oscillations resulting from the homogeneous solution.

There are two important differences when the modulus decay and WIMP
annihilation are included in the evolution. For typical initial values of the
radiation and WIMP, as well as decay rates and branching ratios, the
oscillations in the modulus dominated phase will be over-damped. A second
important effect resulting from decays and annihilations is that these provide
an additional source term which acts to boost the amplitude of the density
perturbations. The enhancement of the amplitude is controlled by the decay and
annihilation rates. This enhancement results in the saturation of the
radiation density perturbation at late times because this growth yields an
additional source for the velocity perturbation resulting in a spatial
dispersion of the radiation fluid. As the radiation velocity perturbation
grows, this slows down the growth of radiation density perturbations. We refer
the reader to Fan:2014zua  for a more detailed analysis.

###  3.2 Evolution through reheating and resulting structure

Once modulus decays significantly affect the energy density (and thus the
expansion rate), this will lead to ‘reheating’ and alter the evolution of the
perturbations. Such effects will become important around an average decay time
$t_{d}\sim H^{-1}_{d}\sim\Gamma_{\phi}^{-1}$. As the modulus decays become
important, the constant scaling from the modulus dominated phase is no longer
valid and we will perform the analysis numerically following the approach in
Fan:2014zua .

The behavior of the radiation perturbations will change near the time of
reheating. Namely, once modulus decay becomes significant it will change the
behavior in which they are enhanced as entering the Hubble radius and then
saturate due to the balance in the source terms.

The modulus decay will typically happen in less than a Hubble time888This does
not have to be the case, but typically is in models resulting from predictions
within particle theory.. This rapid decay results in a termination of the
source terms and the relativistic conversion of modulus particles to radiation
acts to wipe out the growth prior to decay. Thus, the only remnant of the
modulus epoch is an extra suppression in the amplitude of radiation
perturbations (as a consequence of the decay) Fan:2014zua ; Erickcek:2011us ,
which can lead to an enhanced growth of WIMP DM structures on small scales.
The addition of WIMP annihilation as a resulting source for the radiation
density does not change this conclusion. The key effect of adding WIMP
annihilation to the system is to effectively introduce an additional source of
radiation complementary to that provided by modulus decay,999This is
negligible in the branching scenario. which increases the rate of initial
growth and the importance of the saturation by the source terms. This change
also leads to a larger suppression at the time of modulus decay.

For the WIMP perturbations, the density contrast begins to grow linearly with
the scale factor upon horizon entry until reheating. The solution Eq. (3.26)
is an excellent fit to the numerical solution for typical values of $A_{1}$
and $A_{2}$ motivated by particle theory Fan:2014zua . After reheating, the
WIMP annihilation terms become the dominant sources in Eq. (3.5) and, along
with the radiation production, leads to a power loss in the WIMP density
contrast. The annihilation happens in much less than a Hubble time and the
resulting density contrast begins to grow logarithmically, as expected in a
radiation dominated universe Brandenberger:1992dw . In the absence of WIMP
annihilation (i.e. in the branching scenario), this result provides the seeds
for further growth of structure Erickcek:2011us , whereas when annihilation is
important they act to further reduce the strength of the perturbations
following reheating Fan:2014zua .

The WIMP perturbations on scales that enter the Hubble-radius during the EMDE
can experience a significant growth. This growth can lead to formation of
substructure in the form of compact mini-halos or other objects that provide
observations of the non-thermal history Erickcek:2011us ; Fan:2014zua . To
investigate this possibility one needs to take into account cut-off scales
that arise due to kinetic decoupling and free-streaming of DM candidates. We
address these questions next.

Kinetic decoupling and free-streaming of WIMPs

Kinetic equilibrium is maintained between the thermal bath and WIMPs through
elastic and inelastic scattering processes. Once the rate of these scattering
events decreases below the Hubble expansion rate at a temperature $T_{kd}$,
this determines the scale $k_{d}$ of kinetic decoupling.

There are two effects that are important for determining whether the enhanced
growth of perturbations from the matter phase will lead to enhanced growth of
structure on small scales. The kinetic decoupling scale corresponds to the
Hubble scale at the time of decoupling $k^{-1}_{kd}\sim H^{-1}_{kd}$ and at a
temperature $T_{kd}$. However, the process is not instantaneous and for
temperatures below $T_{kd}$, particles can ‘free-stream’ in and out of over-
dense regions in the WIMP distribution. This free-streaming scale is
determined by an integral involving the average velocity following decoupling
until today $t_{0}$ as

$k_{fs}^{-1}=\int_{t_{kd}}^{t_{0}}\frac{\langle v\rangle}{a}dt.$ (3.33)

Moreover, for scales $k^{-1}<k^{-1}_{kd}=1/(a_{kd}H_{kd})$, WIMP particles can
couple to acoustic oscillations of the thermal bath resulting in further
damping on these scales.

To determine which scales maintain the enhancement of the perturbations, it is
useful to introduce a Gaussian cutoff into the perturbations

$\delta_{\chi}=\exp\left[-\frac{k^{2}}{2k_{R}^{2}}\left(\frac{k_{R}}{k_{cut}}\right)^{2}\right]\delta_{\chi}\left(t_{R}\right)$
(3.34)

where

$k_{cut}^{-1}=\mbox{Max}\left(k_{fs}^{-1},k_{kd}^{-1}\right).$ (3.35)

When the modulus decays to WIMPs, the WIMPs will thermalize if their
scattering rate exceeds the Hubble rate. In this case, $T_{kd}<T_{D}^{\phi}$
and $k_{cut}^{-1}=k_{kd}^{-1}$ since $k_{kd}^{-1}>k_{fs}^{-1}$. This leads to
the wash-out of the enhanced growth and is typically the case when non-thermal
DM thermalizes following reheating.

Additionally, WIMPs produced from the modulus decay can have significant
momentum. This can lead to free-streaming of the dark matter for wavelengths
smaller than $k_{fs}^{-1}$. This can also erase the seeded structure resulting
from the enhancement during matter domination and suppress structure
formation. If WIMPs do not thermalize during reheating, then
$T_{kd}>T_{D}^{\phi}$ and the WIMPs will decouple prior to reheating.

There will typically be a continuous momentum spectrum and the average
momentum at reheating $\langle p_{R}\rangle$ is model dependent. However, we
can gain insight by considering when the DM resulting from decay is primarily
non-relativistic $\langle p_{R}\rangle\ll m_{\chi}$. In this case one finds
Fan:2014zua ; Erickcek:2011us

$\displaystyle\frac{k_{R}}{k_{fs}}\approx 2\langle
v_{R}\rangle\left(\sinh^{-1}\sqrt{\frac{\sqrt{2}k_{R}}{k_{eq}}}-\sinh^{-1}\sqrt{a_{eq}}\right),$
(3.36)

whereas for the relativistic case $\langle p_{R}\rangle\gg m_{\chi}$ the
result is

$\displaystyle\frac{k_{R}}{k_{fs}}\approx\frac{a_{NR}}{a_{R}}-1\approx\frac{\langle
p_{R}\rangle}{m_{\chi}},$ (3.37)

where the subscript ‘$eq$’ denotes the time of radiation and matter equality,
and $a_{NR}\equiv\langle p_{R}\rangle a_{R}/m_{\chi}$. From the relativistic
result, one can see that when the decaying modulus is significantly more
massive than the WIMP, free-streaming effects will eliminate the enhanced
growth from the modulus epoch. In the non-relativistic case, however, the
reheat scale may be smaller than the free-streaming scale since

$\displaystyle\frac{k_{rh}}{k_{eq}}=1.2\times
10^{6}\left(\frac{T_{D}^{\phi}}{1\,{\rm
MeV}}\right)\left(\frac{10.75}{g_{*,s}}\right)^{1/3}\left(\frac{g_{*}}{10.75}\right)^{1/2}$
(3.38)

and $\sinh^{-1}x$ behaves as $\log x$ for $x\gg 1$.

##  4 Axion miniclusters and relic density

In this Section, we discuss the formation of axion miniclusters and
expectations for the abundance of axion dark matter in this framework. As we
have previously mentioned, we concentrate on the scenario where PQ breaking
occurs before inflation.

The first-order Boltzmann equations for the modulus and axion, including the
energy exchange terms, are given by Eqs. (3.3), (3.5), (3.7), and (3.9). The
axion perturbation growth is analogous to the modulus case Eq. (3.20):

$\delta_{a}(a,k)=-2\Phi_{0}-\frac{2}{3}\left(\frac{k}{H(T_{\text{osc}}^{a})}\right)^{2}a\Phi_{0}$
(4.1)

where $H(T_{\text{osc}}^{a})$ is the Hubble scale when axion oscillations
begin. We have also defined our scale factor normalization to be
$a(T_{\text{osc}}^{a})\equiv 1$.

In Nelson:2018via , a general set of criteria was defined for the formation of
axion miniclusters which we follow loosely here. Namely, one would expect
miniclusters to form efficiently if a perturbation mode grows to
$\delta_{a}\sim 1$ by the end of the EMD phase. In this work, we also look at
cases where $\delta_{a}=10^{-2}$ and $\delta_{a}=10^{-3}$ which may also be
sufficient for minicluster formation. One might expect the minicluster mass
$M$ at the time of formation can be estimated by Nelson:2018via

$M\sim\frac{4}{3}\pi k^{-3}\,n_{a}(T)\,m_{a}(T_{D}^{\phi})$ (4.2)

for a mode with wavelength $k^{-1}$ entering the horizon at some temperature
$T_{\text{osc}}^{a}\gtrsim T\gtrsim T_{D}^{\phi}$. For modes which are
subhorizon when axion oscillations begin, i.e. for $k>H(T_{\text{osc}}^{a})$,
only miniclusters with very small masses can form despite the fact that
perturbation growth may become large for these modes. Conversely, for
superhorizon modes at the onset of axion oscillations,
$k<H(T_{\text{osc}}^{a})$, sufficient growth of the perturbations may be
difficult. We therefore focus on the modes which enter the horizon close to
$T_{\text{osc}}^{a}$ so that $k\sim H(T_{\text{osc}}^{a})$, which are the
modes mostly likely to grow sufficiently to meet the perturbation growth
criteria while also producing miniclusters with large masses. Taking this into
consideration, Eq. (4.1) can be rewritten as:

$\delta_{a}\simeq-2\Phi_{0}-\frac{2}{3}\left(\frac{H(T_{\text{osc}}^{a})}{H(T_{D}^{\phi})}\right)^{2/3}\Phi_{0}$
(4.3)

where $-\Phi_{0}\sim 10^{-4}$ and we have the conditions
$H(T_{D}^{\phi})\sim\Gamma_{\phi}\propto(T_{D}^{\phi})^{2}$ and
$H(T_{\text{osc}}^{a})\sim m_{a}(T_{\text{osc}}^{a})\propto f_{a}^{-1}$. Thus,
for a given $T_{D}^{\phi}$ we can translate this into a maximal value of the
PQ scale $f_{a}$ which can be expected to form miniclusters with potentially
observable masses.

Figure 3: The maximal value of $f_{a}$ which can be expected to form
miniclusters for a given $T_{D}^{\phi}$. The condition for minicluster
formation is taken to be $\delta_{a}\in\\{1,10^{-2},10^{-3}\\}$ for the solid,
dashed, and dot-dashed curves respectively.

We show in Fig. 3 this bound on the PQ scale for each
$\delta_{a}\in\\{1,10^{-2},10^{-3}\\}$. Here, we note that imposing larger
values of $\delta_{a}$ for a given $f_{a}$ also forces $T_{D}^{\phi}$ to be
closer to the BBN limit. For example, for $f_{a}=10^{9}$ GeV, $T_{D}^{\phi}$
is constrained to be below 33 MeV for $\delta_{a}=1$ while it may be as large
as 240 MeV for $\delta_{a}=10^{-2}$ and 780 MeV for $\delta_{a}=10^{-3}$.
These bounds correspond to $m_{\phi}\lesssim 2.4\times 10^{5}$ GeV,
$m_{\phi}\lesssim 1.2\times 10^{6}$ GeV, and $m_{\phi}\lesssim 2.8\times
10^{6}$ GeV, respectively, for our benchmark effective coupling $c=27.5$. This
behavior can be understood from the fact that achieving larger values of
$\delta_{a}$ require longer periods of matter domination, and therefore lower
$T_{D}^{\phi}$. Conversely, for a given $\delta_{a}$, increasing $f_{a}$ (and
thus lowering the scale of axion oscillations) requires lowering
$T_{D}^{\phi}$ for sufficient overlap. The general trends can be understood as
follows. From Eq. (4.3), it is clear that the perturbation growth scales
parametrically as

$\delta_{a}\,\propto\,\frac{1}{f^{2/3}_{a}(T_{D}^{\phi})^{4/3}}.$ (4.4)

Both qualitative behaviors described above can be understood from the
parametric scaling in Eq. (4.4).

Figure 4: The estimated minicluster mass for modes with $k\sim
H(T_{\text{osc}}^{a})$ at the time of formation, $T_{D}^{\phi}$, for a given
$f_{a}$. Superhorizion modes which have $k\gtrsim H(T_{\text{osc}}^{a})$ may
produce larger miniclusters than those displayed. The condition for
minicluster formation is taken to be $\delta_{a}\in\\{1,10^{-2},10^{-3}\\}$
for the solid, dashed, and dot-dashed curves respectively.

We can also estimate the mass of the miniclusters at the time of formation for
these modes. As we are primarily focused on modes entering the horizon close
to $T_{\text{osc}}^{a}$, Eq. (4.2) can be rewritten as

$M\sim\frac{2}{3}\pi\left(\frac{m_{a}(T_{D}^{\phi})}{m_{a}(T_{\text{osc}}^{a})}\right)\frac{\mathcal{A}_{0}^{2}}{m_{a}(T_{\text{osc}}^{a})}$
(4.5)

where $\mathcal{A}_{0}^{2}$ is given by Eq. (2.29). We show this estimate in
Fig. 4. Here, we adopt $T_{D}^{\phi}$ and corresponding $T_{\text{osc}}^{a}$
values from the bounds on $(T_{D}^{\phi},f_{a})$ pairs for the various
$\delta_{a}$ which are shown in Fig. 3. The differences between the three
curves shown in Fig. 4 all stem from variations of $m(T_{D}^{\phi})$ and
$m(T_{\text{osc}}^{a})$ along these bounds. For values of the PQ scale
$f_{a}\lesssim\mathcal{O}(10^{14}\text{ GeV})$ and
$f_{a}\lesssim\mathcal{O}(10^{10}\text{ GeV})$ for $\delta_{a}=10^{-3}$ and
$\delta_{a}=10^{-2}$, respectively, we have the maximally allowed
$T_{D}^{\phi}\gtrsim\Lambda_{\text{QCD}}$ which causes a reduction in
$m_{a}(T_{D}^{\phi})$ compared to values which have
$T_{D}^{\phi}\lesssim\Lambda_{\text{QCD}}$. Conversely, the larger values of
$T_{D}^{\phi}$ allowed by smaller $\delta_{a}$ require an increase the axion
oscillation temperature, as can be seen from Eq. (2.27). This causes a
reduction in $m_{a}(T_{\text{osc}}^{a})$ which, as can be seen in the figure,
can compete against the previously mentioned reduction in
$m_{a}(T_{D}^{\phi})$. We also note that for $\delta_{a}=1$, we do not display
$f_{a}\gtrsim 10^{16}$ GeV as the minicluster formation criteria can only be
met for $T_{D}^{\phi}$ in violation of BBN bounds.

Finally, it is worth noting that the estimated minicluster masses in Fig. 4
are not necessarily the maximum minicluster masses that one could predict in
this framework. By loosening our demands on the perturbation growth
$\delta_{a}$, it is possible that much larger superhorizon modes
$k<H(T_{\text{osc}}^{a})$ may meet these demands and thus produce miniclusters
with even larger masses. We leave a more precise study of minicluster masses
with weakened perturbation growth constraints and the study of their evolution
to the present day for future work.

Figure 5:  Axion relic density versus $T_{D}^{\phi}$ where $\theta_{i}=2$. For
our benchmark, these values of $T_{D}^{\phi}$ correspond to
$m_{\phi}\in[2.4\times 10^{4},\,5\times 10^{6}]$ GeV. Purple regions are
consistent with $\delta_{a}=1$, dark blue with $\delta_{a}=10^{-2}$, and light
blue with $\delta_{a}=10^{-3}$. Orange regions correspond to even lower values
of $\delta_{a}$ but still have some overlap with an EMDE. Red regions have no
overlap with an EMDE.  Figure 6:  Axion relic density versus $T_{D}^{\phi}$
where $\theta_{i}=3.113$. For our benchmark, these values of $T_{D}^{\phi}$
correspond to $m_{\phi}\in[2.4\times 10^{4},\,5\times 10^{6}]$ GeV. Purple
regions are consistent with $\delta_{a}=1$, dark blue with
$\delta_{a}=10^{-2}$, and light blue with $\delta_{a}=10^{-3}$. Orange regions
correspond to even lower values of $\delta_{a}$ but still have some overlap
with an EMDE. Red regions have no overlap with an EMDE.

Now that we have an understanding of the regions of parameter space where
axion miniclusters may be expected to form, we can look at the expected relic
abundance of axion dark matter within these regions. In Fig. 5, we display the
axion relic density computed from numerical solutions to the Boltzmann
equations as a function of $T_{D}^{\phi}$ for an untuned $\theta_{i}=2$.
Purple regions are consistent with the minicluster formation criterion
$\delta_{a}=1$, dark blue with $\delta_{a}=10^{-2}$, and light blue with
$\delta_{a}=10^{-3}$. Orange regions correspond to even lower values of
$\delta_{a}$ but still have some overlap with an EMDE. Red regions, however,
have $T_{\text{osc}}^{a}<T_{D}^{\phi}$ and thus represent axions which only
begin oscillating after the EMDE phase has ended. From the figure, we also see
the red regions rapidly approach constant values; this is as expected since
for $T_{D}^{\phi}>T_{\text{osc}}^{a}$, the axion is still frozen until any
entropy release from modulus decay has ceased. However, any overlap with an
EMDE then necessarily implies overlap with the era of entropy release,
resulting in a reduction of the axion yield $Y_{a}$ and thus a reduction in
its relic abundance.

It is immediately evident from Fig. 5 that, at least for EMDEs consistent with
our demand of minicluster formation, saturating the total DM relic density
with axions is only viable for $f_{a}\gtrsim\mathcal{O}(10^{13})$ GeV. For
lower values of $f_{a}$, the horizontal dashed line displaying $\Omega
h^{2}=0.12$ can only intersect with red regions. Additionally, if we impose
the demand for both minicluster formation and a DM relic density which is
saturated by axions, we are forced to have
$T_{D}^{\phi}\,\lesssim\mathcal{O}(100)$ MeV. As expected, the regions
compatible with minicluster formation become smaller the larger $\delta_{a}$
is. For $\delta_{a}=1$, for example, we are forced into the purple regions
which only intersect the relic density line when $\phi$ decays very close to
BBN and $f_{a}\sim\mathcal{O}(10^{16})$ GeV. Additionally, if $T_{D}^{\phi}$
is further required to be above $\sim 4$ MeV for compatibility with BBN as
suggested in Kawasaki:2000en ; Hasegawa:2019jsa , our results become even more
constrained so that only $f_{a}\lesssim\mathcal{O}(10^{14})$ GeV is allowed
and the total relic density cannot be satisfied with the choice $\theta_{i}=2$
for the misalignment angle.

The general trends can once again be understood from the parametric scaling in
Eq. (4.4). For a given value of $f_{a}$, increasing $T_{D}^{\phi}$ means
shorter periods of matter domination and hence smaller values of $\delta_{a}$;
thus, purple regions transition to dark blue, then light blue, and finally
orange. On the other hand, for a given $T_{D}^{\phi}$ increasing $f_{a}$ also
lowers $\delta_{a}$, causing the viable regions for minicluster formation to
shrink as $f_{a}$ increases.

For values of the misalignment angle which approach the limiting value,
$\theta_{i}\rightarrow\pi$, the axion relic abundance can be enhanced by a
potentially large amount due to anharmonic effects Visinelli:2009zm . We
display the axion relic density for a tuned $\theta_{i}=3.113$ in Fig. 6,
again computed from numerical solutions to the Boltzmann equations. The color
coding is identical to Fig. 5. Generally, the trends discussed in the case of
$\theta_{i}=2$ continue to hold, with the contours of constant $f_{a}$ pushed
upwards to higher relic densities. It is clear that, even with a mild amount
of tuning, demanding both minicluster formation and purely axion DM is only
viable for low $T_{D}^{\phi}\lesssim\mathcal{O}(200)$ MeV and
$f_{a}\sim\mathcal{O}(10^{13})$ GeV. Meeting both demands for lower values of
the PQ scale $f_{a}$ thus clearly requires a large degree of tuning for the
misalignment angle.

##  5 Relic density contributions from WIMPs

In light of the results in Section 4, we might expect that if axion
miniclusters do form, there should be additional contributions from other
species to the total amount of dark matter. The parameter space with
$\delta_{a}\sim\mathcal{O}(10^{-3}-1)$ for values of
$f_{a}\lesssim\mathcal{O}(10^{14})$ GeV typically has axions contributing as a
sub-dominant component of dark matter, as seen in both Fig. 5 for an untuned
case and in Fig. 6 for a mildly tuned case. Demanding values of $\delta_{a}$
at the larger end of this range forces one into larger values of $f_{a}$ if
axions are to constitute a significant percentage of the total dark matter.
Indeed, if one requires $\delta_{a}\sim 1$, one ultimately runs out of
parameter space even for large $f_{a}$, as the BBN limit on $T_{D}^{\phi}$ is
breached.

Since it appears that axion-only dark matter is unlikely to saturate the
observed dark matter density while also producing miniclusters, we are led to
consider a multi-component dark matter scenario. An appreciable amount - if
not a majority - of the dark matter might then be composed of some other dark
matter particle, which allows for the formation of axion miniclusters while
easing the requirement of extreme misalignment angle tuning. In this work, we
study the case where the second component is composed of WIMPs.

Generally, one would expect that the WIMP annihilation cross section can be
well-approximated by $s$-wave and $p$-wave contributions:

$\langle\sigma v\rangle_{\chi}\sim s_{0}+p_{0}\frac{T}{m_{\chi}}.$ (5.1)

From the WIMP Boltzmann equation given in Eq. (2.14), we can see that while
WIMPs are relativistic, the annihilation cross section primarily serves to
maintain equilibrium. However, once $T\lesssim m_{\chi}$ and the equilibrium
density $\overline{\rho}_{\chi}$ becomes Boltzmann-suppressed, the precise
value of $\langle\sigma v\rangle_{\chi}$ plays a large role in determining the
final WIMP relic density. Thus, the $p$-wave contributions are only very
important for our considerations if $p_{0}\gtrsim s_{0}$. Here we focus on the
case of a purely $s$-wave annihilation cross section so that $\langle\sigma
v\rangle_{\chi}=s_{0}$ and leave the analysis of the $p$-wave contributions
for future work.

WIMP annihilations can produce detectable gamma-ray signals. Such signals have
been the focus of many experiments such as HESS HESS:2016mib , MAGIC
MAGIC:2016xys , and Fermi-LAT Fermi-LAT:2015att ; Fermi-LAT:2016uux , which
have constrained the $T\rightarrow 0$ limit of $\langle\sigma v\rangle_{\chi}$
for several different production channels. Here we take a model-independent
approach and adopt the limits from the $\tau^{+}\tau^{-}$ channel shown in
Fermi-LAT:2015att . This particular channel tends to be an upper bound on the
WIMP annihilation cross section. Other channels (if relevant for a given
model) will then be even more highly constrained. As we will see, even under
the most optimistic constraints it is extremely difficult to produce
reasonable amounts of WIMP dark matter in the regime where axion miniclusters
may be expected to form.

Figure 7: Upper limits on the annihilation of WIMP DM to $\tau^{+}\tau^{-}$
taken from Fermi-LAT Fermi-LAT:2015att . Here, we rescale the constraints
assuming that the observed DM density is composed of 100% WIMPs (blue curve),
10% WIMPs (orange curve), and 1% WIMPs (green curve).

Firstly, the Fermi-LAT upper bounds on the annihilation cross section to
$\tau^{+}\tau^{-}$ needs to be reinterpreted in the context of a multi-
component dark matter scenario. The constraints become weaker if WIMPs make up
a smaller percentage of the total dark matter density; specifically, the WIMP
constraint $\langle\sigma v\rangle_{\chi}^{\text{multi}}$ in a multi-component
scenario can be related to the constraint $\langle\sigma
v\rangle_{\chi}^{\text{single}}$ in a WIMP-only scenario by Baer:2016ucr

$\langle\sigma v\rangle_{\chi}^{\text{multi}}=\xi^{2}\langle\sigma
v\rangle_{\chi}^{\text{single}}$ (5.2)

where $\xi\leq 1$ is the percentage of total dark matter being composed of
WIMPs. In Fig. 7, we show the upper limits on the annihilation of WIMPs from
Fermi-LAT, appropriately rescaled assuming that the observed DM density is
composed of 100% WIMPs (blue curve), 10% WIMPs (orange curve), and 1% WIMPs
(green curve).

Figure 8:  Annihilation Scenario: Normalized WIMP relic density
$\Omega_{\chi}h^{2}$ versus WIMP mass $m_{\chi}$ for $m_{\phi}=4.1\times
10^{5}$ GeV ($T^{\phi}_{D}=51$ MeV). We display curves for values of
$\xi\in\\{0.001,\,0.01,\,0.1,\,0.5,\,1\\}$ where $\xi$ is the percentage of
total dark matter composed of WIMPs. We fix the annihilation cross section
$\langle\sigma v\rangle_{\chi}^{\text{multi}}$ to the maximally allowed value
by Fermi-LAT constraints, $\xi^{2}\langle\sigma
v\rangle_{\chi}^{\text{single}}$, as in Fig. 7. We have also fixed
$B_{\phi\rightarrow\chi}=5\%$.  Figure 9:  Branching Scenario: We display
contours of the branching of the modulus to WIMPs $B_{\phi\rightarrow\chi}$
for which the observed DM relic density is obtained, assuming Fermi-LAT
constraints with $\xi=1$. Here, the solid boundaries of contours show where
the overall DM relic density is saturated exclusively by WIMPs, while shaded
regions have $\Omega_{\chi}h^{2}<0.12$. For $T_{D}^{\phi}\gtrsim 0.2$ GeV, the
scenario transitions to the annihilation scenario for the values of
$B_{\phi\rightarrow\chi}$ shown. The black vertical lines indicate the highest
$T^{\phi}_{D}$ for which axion minicluster formation is expected to occur for
$\delta_{a}\in\\{1,\,10^{-2},\,10^{-3}\\}$.

We now discuss the confluence of the axion perturbation condition, relic
density, and Fermi-LAT upper limits on the WIMP annihilation cross section.
From Sec. 4 we see that the axion perturbation condition prefers low
$T^{\phi}_{D}\,\sim\,\mathcal{O}(1-10)$ MeV. We first understand the
implications of this for WIMPs in the annihilation scenario, where the relic
density is given by Eq. (2.24). For low $T^{\phi}_{D}$ and
$T_{f}\,\sim\,m_{\chi}/20\,\sim\,5$ GeV (for a 100 GeV WIMP), we have

$T_{f}/T^{\phi}_{D}\,\sim\,\mathcal{O}(100)$ (5.3)

which implies that the annihilation scenario will generally oversaturate the
relic density if the condition for axion miniclusters is satisfied. To
compensate for the large factor in Eq. (5.3) we consider WIMP candidates with
highly enhanced annihilation cross sections; however, that option is subject
to upper limits from Fermi-LAT.

This situation is portrayed in Fig. 8. We take $B_{\phi\rightarrow\chi}=5\%$
so that we are indeed in the annihilation scenario. We also set
$m_{\phi}=4.1\times 10^{5}$ GeV so that $T^{\phi}_{D}\simeq 51$ MeV, roughly
corresponding to the upper limit allowed by $\delta_{a}=1$. The vertical axis
is the normalized WIMP contribution to the DM relic density, which should be -
at most \- unity (shown by the black dashed line) for consistency. The blue
curve depicts the normalized WIMP contribution to the relic density
corresponding to the Fermi-LAT upper limits on the annihilation cross section
when WIMPs comprise 100% of DM. Clearly, WIMPs oversaturate the relic density
in this case, unless they are extremely heavy and above the reach of Fermi-
LAT. Reducing the WIMP contribution (reducing $\xi$) mitigates the tension
somewhat, but is generally unable to overcome the effect of the large value of
$T_{f}/T^{\phi}_{D}$. We therefore see that if the annihilation cross section
must be raised significantly to produce a sensible quantity of WIMPs, we are
forced to consider scenarios where WIMPs make up a small percentage of the
total dark matter. We are therefore again plagued by the tuning problems of
our axion-only DM scenario of the previous section unless yet another DM
candidate is considered.

We now move on to the branching scenario for WIMPs. In this case, the WIMP
yield is given by Eq. (2.26). We first describe the parametric scaling of the
DM yield and modulus branching, and Fig. 9 gives our quantitative results.
Firstly, we note that the modulus yield scales parametrically as

$Y_{\phi}\,\sim\,T^{\phi}_{D}/m_{\phi}\,\propto\,(T^{\phi}_{D})^{1/3}\,\,.$
(5.4)

Comparing to dark matter observations implies a WIMP yield of

$Y_{\chi}\simeq 5\times 10^{-10}\left(\frac{1\text{ GeV}}{m_{\chi}}\right)$
(5.5)

and so the branching ratio is

$B_{\phi\rightarrow\chi}\sim\mathcal{O}\left(10^{-4}\right)\times\left(\frac{10\text{
GeV}}{m_{\chi}}\right)\left(\frac{m_{\phi}}{10^{5}\text{
GeV}}\right)\left(\frac{0.05\text{ GeV}}{T_{D}^{\phi}}\right).$ (5.6)

Thus, for $T^{\phi}_{D}\sim\mathcal{O}(1-50)$ MeV required for axion
minicluster formation the branching of the modulus to WIMPs falls in the range
$B_{\phi\rightarrow\chi}\sim 10^{-4}-10^{-7}$ depending on the precise
magnitude of $m_{\chi}$. This is shown in Fig. 9. The larger branching ratios
are concentrated on the top right corner, where the annihilation scenario is
operational. For the smaller branchings displayed $B_{\phi\rightarrow\chi}\sim
10^{-4}-10^{-6}$, the annihilation scenario is also operational beginning at
$T_{D}^{\phi}\gtrsim 0.2$ GeV for small WIMP masses. This results from
$m_{\chi}$ being sufficiently low so that the late WIMP freeze-out temperature
(close to $T_{D}^{\phi}$) results in $\overline{\rho}_{\chi}$ being large
enough to produce WIMPs through thermal interactions. For larger $m_{\chi}$,
the WIMP freeze-out temperature is then larger than $T_{D}^{\phi}$ so these
reactions do not occur and the branching scenario takes over.

In Fig. 9, the vertical dashed lines show the maximal value of $T^{\phi}_{D}$
that is compatible with axion minicluster formation, as obtained from Fig. 5.
We see that for the $\delta_{a}=1$ bound, the total DM relic density can
indeed be saturated for the displayed branching ratios
$B_{\phi\rightarrow\chi}\sim 10^{-4}-10^{-6}$ for most of the range of
$m_{\chi}$ probed by Fermi-LAT. However as $\delta_{a}$ decreases the
parameter space becomes more restricted, requiring heavier WIMPs with
$m_{\chi}\gtrsim 1-2$ TeV to have $B_{\phi\rightarrow\chi}\lesssim 10^{-7}$.
Many currently existing models involving moduli and WIMPs predict
$B_{\phi\rightarrow\chi}\sim\mathcal{O}(0.01-0.1)$, and are thus in tension
with these results - although this certainly does not preclude the possibility
of a compatible WIMP model with such a low $B_{\phi\rightarrow\chi}$. We
reiterate that here we take an agnostic approach as to the model-building
details and simply regard the branching ratio as a free parameter.

Based on these results, it appears that the branching scenario - which is
capable of saturating the total DM relic density - may be preferable in
regions of parameter space that allow for axion minicluster formation within a
multicomponent WIMP and axion dark matter framework. It is interesting to note
that, as discussed in Sec. 3, the branching scenario may allow for the
sustained growth of WIMP DM substructures while the annihilation scenario
typically has any growth washed out by rapid thermalization. Thus, unless
these effects are washed out by free-streaming (as may be likely for small
WIMP masses), it appears plausible that WIMP substructures are likely form in
tandem with axion miniclusters within this framework. We leave a thorough
analysis of these substructures for future work.

##  6 Conclusions

In this paper, we have studied the question of DM substructure formation in
the presence of an EMDE, with a focus on axion miniclusters. We have
especially dealt with the case of PQ breaking before inflation, where EMDEs
can play a critical role in ensuring minicluster formation. Our results are as
follows.

Firstly, the axion oscillation temperature must be larger than the modulus
decay temperature: $T^{a}_{\rm osc}\,>\,T^{\phi}_{D}$. Secondly, the axion
minicluster formation condition Eq. (4.3) requires
$H(T_{\text{osc}}^{a})\,\gg\,H(T_{D}^{\phi})$. However, since
$H(T_{D}^{\phi})\sim\Gamma_{\phi}\propto(T_{D}^{\phi})^{2}$ and
$H(T_{\text{osc}}^{a})\sim m_{a}(T_{\text{osc}}^{a})\propto f_{a}^{-1}$, these
conditions effectively force $T^{\phi}_{D}$ to be near BBN, and especially for
cases when the axion constitutes a significant part of the relic density. The
effect of demanding axion minicluster on WIMPs then is that they are forced
into the branching scenario, in order to be compatible with upper limits
coming from Fermi-LAT. If WIMPs comprise a significant percentage of the total
dark matter so the constraints on the PQ sector are eased, a branching ratio
$B_{\phi\rightarrow\chi}\lesssim\mathcal{O}(10^{-4})$ is then required so that
the total relic density does not become oversaturated. Furthermore, because
the branching scenario is required in this regime, our results suggest that
models which can meet the demands on the PQ sector and on
$B_{\phi\rightarrow\chi}$ for minicluster formation may also predict WIMP DM
substructures - especially for heavier WIMPs - unless free-streaming effects
become significant.

Despite the effective decoupling of WIMPs and axions, the level of constraint
on one sector by imposing demands on the other within our EMDE framework is
quite remarkable. These conclusions also pave several avenues for future work.
In particular, a detailed study of free-streaming effects for WIMPs produced
through the branching scenario would reveal to what degree WIMP substructure
formation is indeed correlated with axion substructure formation.
Additionally, the inclusion of sizeable $p$-wave contributions to the WIMP
cross section may change these predictions as 1.) the WIMP relic density may
be decreased by a large degree in the annihilation scenario while the
$T\rightarrow 0$ limit on $\langle\sigma v\rangle_{\chi}$ remains compatible
with Fermi-LAT bounds, and 2.) this may result in a larger degree of WIMP
thermalization in the branching scenario, which likely erases any tentative
WIMP substructure growth. Finally, it is interesting to question whether other
dark matter candidates such as Primordial Black Holes could share such a high
degree of constraint in the context of a similar EMDE multicomponent DM model.

## Acknowledgements

We thank Amy Burks, Adrienne Erickeck, and Nemanja Kaloper for useful
conversations. S.W. thanks KITP Santa Barbara and the Simons Center for
hospitality. This research was supported in part by DOE grant DE-
FG02-85ER40237. K.S. thanks the Simons Center for hospitality and is supported
in part by DOE grant DE-SC0009956. The computing for this project was
performed at the OU Supercomputing Center for Education & Research (OSCER) at
the University of Oklahoma (OU).

## References

  * (1) S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * (2) R.D. Peccei and H.R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * (3) F. Wilczek, _Problem of Strong p and t Invariance in the Presence of Instantons_ , _Phys. Rev. Lett._ 40 (1978) 279.
  * (4) J.E. Kim, _Weak Interaction Singlet and Strong CP Invariance_ , _Phys. Rev. Lett._ 43 (1979) 103.
  * (5) M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, _Can Confinement Ensure Natural CP Invariance of Strong Interactions?_ , _Nucl. Phys. B_ 166 (1980) 493.
  * (6) A.R. Zhitnitsky, _On Possible Suppression of the Axion Hadron Interactions. (In Russian)_ , _Sov. J. Nucl. Phys._ 31 (1980) 260.
  * (7) M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137.
  * (8) L.J. Hall, J.D. Lykken and S. Weinberg, _Supergravity as the Messenger of Supersymmetry Breaking_ , _Phys. Rev. D_ 27 (1983) 2359.
  * (9) B. de Carlos, J.A. Casas, F. Quevedo and E. Roulet, _Model independent properties and cosmological implications of the dilaton and moduli sectors of 4-d strings_ , _Phys. Lett. B_ 318 (1993) 447 [hep-ph/9308325].
  * (10) G.D. Coughlan, W. Fischler, E.W. Kolb, S. Raby and G.G. Ross, _Cosmological Problems for the Polonyi Potential_ , _Phys. Lett. B_ 131 (1983) 59.
  * (11) T. Banks, D.B. Kaplan and A.E. Nelson, _Cosmological implications of dynamical supersymmetry breaking_ , _Phys. Rev. D_ 49 (1994) 779 [hep-ph/9308292].
  * (12) R. Allahverdi et al., _The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe_ , 2006.16182.
  * (13) G. Kane, K. Sinha and S. Watson, _Cosmological Moduli and the Post-Inflationary Universe: A Critical Review_ , _Int. J. Mod. Phys. D_ 24 (2015) 1530022 [1502.07746].
  * (14) G.F. Giudice, E.W. Kolb and A. Riotto, _Largest temperature of the radiation era and its cosmological implications_ , _Phys. Rev. D_ 64 (2001) 023508 [hep-ph/0005123].
  * (15) G.B. Gelmini and P. Gondolo, _Neutralino with the right cold dark matter abundance in (almost) any supersymmetric model_ , _Phys. Rev. D_ 74 (2006) 023510 [hep-ph/0602230].
  * (16) R. Allahverdi, M. Cicoli, B. Dutta and K. Sinha, _Nonthermal dark matter in string compactifications_ , _Phys. Rev. D_ 88 (2013) 095015 [1307.5086].
  * (17) R. Allahverdi, M. Cicoli, B. Dutta and K. Sinha, _Correlation between Dark Matter and Dark Radiation in String Compactifications_ , _JCAP_ 10 (2014) 002 [1401.4364].
  * (18) M. Cicoli, K. Sinha and R. Wiley Deal, _The dark universe after reheating in string inflation_ , _JHEP_ 12 (2022) 068 [2208.01017].
  * (19) B. Dutta, L. Leblond and K. Sinha, _Mirage in the Sky: Non-thermal Dark Matter, Gravitino Problem, and Cosmic Ray Anomalies_ , _Phys. Rev. D_ 80 (2009) 035014 [0904.3773].
  * (20) B.S. Acharya, P. Kumar, K. Bobkov, G. Kane, J. Shao and S. Watson, _Non-thermal Dark Matter and the Moduli Problem in String Frameworks_ , _JHEP_ 06 (2008) 064 [0804.0863].
  * (21) B.S. Acharya, G. Kane, S. Watson and P. Kumar, _A Non-thermal WIMP Miracle_ , _Phys. Rev. D_ 80 (2009) 083529 [0908.2430].
  * (22) C. Blanco, M.S. Delos, A.L. Erickcek and D. Hooper, _Annihilation Signatures of Hidden Sector Dark Matter Within Early-Forming Microhalos_ , _Phys. Rev. D_ 100 (2019) 103010 [1906.00010].
  * (23) R. Easther, R. Galvez, O. Ozsoy and S. Watson, _Supersymmetry, Nonthermal Dark Matter and Precision Cosmology_ , _Phys. Rev. D_ 89 (2014) 023522 [1307.2453].
  * (24) G. Kane, J. Shao, S. Watson and H.-B. Yu, _The Baryon-Dark Matter Ratio Via Moduli Decay After Affleck-Dine Baryogenesis_ , _JCAP_ 11 (2011) 012 [1108.5178].
  * (25) R. Allahverdi, B. Dutta and K. Sinha, _Cladogenesis: Baryon-Dark Matter Coincidence from Branchings in Moduli Decay_ , _Phys. Rev. D_ 83 (2011) 083502 [1011.1286].
  * (26) R. Allahverdi, B. Dutta and K. Sinha, _Baryogenesis and Late-Decaying Moduli_ , _Phys. Rev. D_ 82 (2010) 035004 [1005.2804].
  * (27) R. Allahverdi, B. Dutta and K. Sinha, _Successful Supersymmetric Dark Matter with Thermal Over/Under-Abundance from Late Decay of a Visible Sector Scalar_ , _Phys. Rev. D_ 87 (2013) 075024 [1212.6948].
  * (28) C. Miller, A.L. Erickcek and R. Murgia, _Constraining nonthermal dark matter’s impact on the matter power spectrum_ , _Phys. Rev. D_ 100 (2019) 123520 [1908.10369].
  * (29) M.S. Delos, A.L. Erickcek, A.P. Bailey and M.A. Alvarez, _Density profiles of ultracompact minihalos: Implications for constraining the primordial power spectrum_ , _Phys. Rev. D_ 98 (2018) 063527 [1806.07389].
  * (30) M.S. Delos, A.L. Erickcek, A.P. Bailey and M.A. Alvarez, _Are ultracompact minihalos really ultracompact?_ , _Phys. Rev. D_ 97 (2018) 041303 [1712.05421].
  * (31) J. Fan, O. Özsoy and S. Watson, _Nonthermal histories and implications for structure formation_ , _Phys. Rev. D_ 90 (2014) 043536 [1405.7373].
  * (32) H. Baer, V. Barger, S. Salam, D. Sengupta and K. Sinha, _Status of weak scale supersymmetry after LHC Run 2 and ton-scale noble liquid WIMP searches_ , _Eur. Phys. J. ST_ 229 (2020) 3085 [2002.03013].
  * (33) H. Baer, V. Barger, S. Salam, H. Serce and K. Sinha, _LHC SUSY and WIMP dark matter searches confront the string theory landscape_ , _JHEP_ 04 (2019) 043 [1901.11060].
  * (34) K.J. Bae, H. Baer and E.J. Chun, _Mixed axion/neutralino dark matter in the SUSY DFSZ axion model_ , _JCAP_ 12 (2013) 028 [1309.5365].
  * (35) H. Baer, A. Lessa and W. Sreethawong, _Coupled Boltzmann calculation of mixed axion/neutralino cold dark matter production in the early universe_ , _JCAP_ 01 (2012) 036 [1110.2491].
  * (36) K.J. Bae, H. Baer, A. Lessa and H. Serce, _Coupled Boltzmann computation of mixed axion neutralino dark matter in the SUSY DFSZ axion model_ , _JCAP_ 10 (2014) 082 [1406.4138].
  * (37) A.E. Nelson and H. Xiao, _Axion Cosmology with Early Matter Domination_ , _Phys. Rev. D_ 98 (2018) 063516 [1807.07176].
  * (38) N. Blinov, M.J. Dolan and P. Draper, _Imprints of the Early Universe on Axion Dark Matter Substructure_ , _Phys. Rev. D_ 101 (2020) 035002 [1911.07853].
  * (39) L. Visinelli and J. Redondo, _Axion Miniclusters in Modified Cosmological Histories_ , _Phys. Rev. D_ 101 (2020) 023008 [1808.01879].
  * (40) M.P. Hertzberg, M. Tegmark and F. Wilczek, _Axion Cosmology and the Energy Scale of Inflation_ , _Phys. Rev. D_ 78 (2008) 083507 [0807.1726].
  * (41) D.J. Gross, R.D. Pisarski and L.G. Yaffe, _QCD and Instantons at Finite Temperature_ , _Rev. Mod. Phys._ 53 (1981) 43.
  * (42) M. Dine, L. Randall and S.D. Thomas, _Baryogenesis from flat directions of the supersymmetric standard model_ , _Nucl. Phys. B_ 458 (1996) 291 [hep-ph/9507453].
  * (43) M. Dine, L. Randall and S.D. Thomas, _Supersymmetry breaking in the early universe_ , _Phys. Rev. Lett._ 75 (1995) 398 [hep-ph/9503303].
  * (44) M. Cicoli, K. Dutta, A. Maharana and F. Quevedo, _Moduli Vacuum Misalignment and Precise Predictions in String Inflation_ , _JCAP_ 08 (2016) 006 [1604.08512].
  * (45) V. Poulin, T.L. Smith, D. Grin, T. Karwal and M. Kamionkowski, _Cosmological implications of ultralight axionlike fields_ , _Phys. Rev. D_ 98 (2018) 083525 [1806.10608].
  * (46) L. Visinelli and P. Gondolo, _Axion cold dark matter in non-standard cosmologies_ , _Phys. Rev. D_ 81 (2010) 063508 [0912.0015].
  * (47) R. Hlozek, D.J.E. Marsh and D. Grin, _Using the Full Power of the Cosmic Microwave Background to Probe Axion Dark Matter_ , _Mon. Not. Roy. Astron. Soc._ 476 (2018) 3063 [1708.05681].
  * (48) E.W. Kolb and I.I. Tkachev, _Axion miniclusters and Bose stars_ , _Phys. Rev. Lett._ 71 (1993) 3051 [hep-ph/9303313].
  * (49) C.J. Hogan and M.J. Rees, _AXION MINICLUSTERS_ , _Phys. Lett. B_ 205 (1988) 228.
  * (50) M.S. Turner, F. Wilczek and A. Zee, _Formation of Structure in an Axion Dominated Universe_ , _Phys. Lett. B_ 125 (1983) 35.
  * (51) L. Visinelli and P. Gondolo, _Dark Matter Axions Revisited_ , _Phys. Rev. D_ 80 (2009) 035024 [0903.4377].
  * (52) Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * (53) M. Kawasaki, K. Kohri and N. Sugiyama, _Cosmological constraints on late time entropy production_ , _Phys. Rev. Lett._ 82 (1999) 4168 [astro-ph/9811437].
  * (54) M. Kawasaki, K. Kohri and N. Sugiyama, _MeV scale reheating temperature and thermalization of neutrino background_ , _Phys. Rev. D_ 62 (2000) 023506 [astro-ph/0002127].
  * (55) T. Hasegawa, N. Hiroshima, K. Kohri, R.S.L. Hansen, T. Tram and S. Hannestad, _MeV-scale reheating temperature and thermalization of oscillating neutrinos by radiative and hadronic decays of massive particles_ , _JCAP_ 12 (2019) 012 [1908.10189].
  * (56) K.J. Bae, H. Baer, V. Barger and R.W. Deal, _The cosmological moduli problem and naturalness_ , _JHEP_ 02 (2022) 138 [2201.06633].
  * (57) R. Wiley Deal, _The dark universe: the interplay of cosmological moduli, axions, and the MSSM_ , Ph.D. thesis, Oklahoma U., 5, 2023.
  * (58) E.W. Kolb and M.S. Turner, _The Early Universe_ , vol. 69 (1990), 10.1201/9780429492860.
  * (59) H. Baer, V. Barger and R. Wiley Deal, _Dark matter and dark radiation from the early universe with a modulus coupled to the PQMSSM_ , 2301.12546.
  * (60) P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau et al., _SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python_ , _Nature Methods_ 17 (2020) 261.
  * (61) A.L. Erickcek and K. Sigurdson, _Reheating Effects in the Matter Power Spectrum and Implications for Substructure_ , _Phys. Rev. D_ 84 (2011) 083503 [1106.0536].
  * (62) R.H. Brandenberger, H. Feldman, V.F. Mukhanov and T. Prokopec, _Gauge invariant cosmological perturbations: Theory and applications_ , .
  * (63) H.E.S.S. collaboration, _Search for dark matter annihilations towards the inner Galactic halo from 10 years of observations with H.E.S.S_ , _Phys. Rev. Lett._ 117 (2016) 111301 [1607.08142].
  * (64) MAGIC, Fermi-LAT collaboration, _Limits to Dark Matter Annihilation Cross-Section from a Combined Analysis of MAGIC and Fermi-LAT Observations of Dwarf Satellite Galaxies_ , _JCAP_ 02 (2016) 039 [1601.06590].
  * (65) Fermi-LAT collaboration, _Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data_ , _Phys. Rev. Lett._ 115 (2015) 231301 [1503.02641].
  * (66) Fermi-LAT, DES collaboration, _Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-LAT_ , _Astrophys. J._ 834 (2017) 110 [1611.03184].
  * (67) H. Baer, V. Barger and H. Serce, _SUSY under siege from direct and indirect WIMP detection experiments_ , _Phys. Rev. D_ 94 (2016) 115019 [1609.06735].