Multi-component Dark Matter and Small Scale Structure Formation

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})$.


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 [1][2][3][4][5][6][7].
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 [8][9][10][11], 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 [12,13]; for earlier studies, we refer to [14,15]) and their connections to the physics of moduli stabilization [16][17][18][19][20][21].
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 [22,23], connections with baryogenesis [24][25][26][27] and structure formation [28][29][30][31] 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 substructure 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 [32][33][34][35][36].
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) 1 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 [37][38][39].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 substructure formation prefers low modulus reheat temperatures T ϕ D ≲ O(50) MeV, which can only accommodate pure axion DM for decay constants f a ≳ 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 ∼ 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 ≲ O(10 −4 ), corresponding to the "branching scenario" [12,13].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.

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.

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 ξ is given by where f a is the vacuum expectation value of the PQ field.The expectation value of the ξ can be written as ⟨ξ(x)⟩ = f a e iφ(x)/fa where the phase field φ(x) is identified as the axion field.Here, the axion field φ(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 form2 [40,41] where b = 0.018 and Λ QCD ∼ 200 MeV.The axion potential after PQ symmetry breaking can then be written approximately as 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 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 (∼ 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 φ, behaves like a time-dependent friction term which does not allow the axion field to oscillate until H ≃ m a .This gives the condition for critical time t osc at which the axion oscillations begin: (2.5) We can thus take the axion field to be frozen for t < t osc , which produces an equation of state w a = −1.The axion field then oscillates very rapidly for t > t osc such that its average equation of state in this regime is ⟨w a ⟩ = 0.In other words, the axion field behaves similarly to a cosmological constant before t osc with its number density n a = constant 3 , and it behaves like dark matter after t osc where its number density dilutes as n a ∼ a −3 where we use a to denote the scale factor.Thus, we can estimate the number density of axions as where θ i = φ 0 /f a ∈ [0, π) is the initial misalignment angle, which is random and is typically assumed to be O(1).Now, the axion energy density at time t is simply given by ρ 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 [42,43] and frozen in place by Hubble friction.Once the universe has expanded sufficiently (i.e. when H ∼ m ϕ ), 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 [42][43][44].Additionally, moduli are typically assumed to be very massive with m ϕ ≳ O(10 − 50 TeV) to avoid conflicts with BBN.

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 [35,39,40,[45][46][47].Denoting the Hubble scale during inflation by H I , the two scenarios are f a < H I /2π (axion absent during inflation) and f a > H I /2π (axion present during inflation).For the former case, PQ symmetry breaking occurs after inflation and the initial misalignment angle θ i acquires different values within the same Hubble patch.The formation of axion miniclusters in this case has been studied extensively [48][49][50] 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 , θ i , and f a as follows [51,52]: Constraints coming from isocurvature fluctuations in the case of scenarios with decaying scalar fields were worked out in [46].In this case, one has which holds for f a < 10 15 GeV.For higher values of f a ≳ 10 15 GeV, one has GeV . (2.9) For details of other regimes of f a , we refer to [46].Independent of axion physics, the scale of inflation is constrained by Planck [52] to be H I ≲ 2.5 × 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.

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 with the resulting Hubble equation where α ∈ {ϕ, χ, r, a} represents moduli, WIMPs, radiation, and axions, respectively.We work in units of the reduced Planck mass m p = 2.44 × 10 18 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 ) ) subject to the energy constraint Eq. (2.12) and where we have kept the equation of state w α 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. 4Indeed, 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.

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 (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 as 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 ≲ O(1 MeV), decay occurs during or after Big Bang Nucleosynthesis (BBN) [53][54][55] 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 Γ ϕ ≳ 5 × 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 [55], 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 .From imposing ρ r ∼ ρ ϕ and assuming a radiation-dominated universe before modulus domination, the temperature of modulus-radiation equality can be shown to be 5 [56] where ϕ 0 is the initial amplitude of the modulus.As the universe is matter-dominated between T ϕ e and T ϕ D , it is a simple matter to relate the Hubble scales at both eras.Upon noting that entropy is conserved between T ϕ e and -by definition -the entropy injection temperature T ϕ S , while entropy is not conserved between T ϕ S and T ϕ D , we can arrive at the following estimate for the entropy injection temperature: 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 freezeout scale and assuming the WIMPs behave as cold dark matter in this regime (i.e.w χ = 0 and ρ χ ≃ m χ n χ ), the Boltzmann equation for WIMPs can be rewritten in terms of the abundance yield Taking the limit of constant ⟨σv⟩ χ , this equation can be straightforwardly integrated from the scale of modulus decay T ϕ D to the present time to give the result Assuming no significant source of entropy production after the modulus has decayed, the WIMP relic density can be estimated by [58] where s 0 ∼ 2970 cm −3 and ρ c /h 2 ∼ 1.05 × 10 −5 GeV • cm −3 .This can also be used to relate the produced WIMP relic density Ω χ h 2 to the expected WIMP relic density produced in a standard thermal scenario This well-known result is typically coined the annihilation scenario in the literature [16][17][18].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 where we again rewrite the Boltzmann equation in terms of the abundance yield Y χ .Integrating this equation from the scale of modulus decay to the present time then gives the result Once again, the WIMP relic density can be estimated with Eq. (2.23).This result is typically denoted as the branching scenario [16][17][18].
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 < T a osc , so that axion miniclusters might form if PQ symmetry breaking occurs at or above the inflationary scale [37].This has the crucial property that entropy dilution from the modulus will reduce the late-time relic density of the misalignment axions [59].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 [58] where we assume a transition to a radiation-dominated universe at T ϕ D .It is then straightforward to compute T a osc from Eq. (2.2) and Eq.(2.5), which we find to be roughly T a osc ∼ O(0.1 − 1 GeV) depending on the value of T ϕ D .Having found the axion oscillation temperature, we can then estimate the axion number density when oscillations begin by where, accounting for anharmonic effects, the initial oscillation amplitude A 0 is given by [51] (2.29) in terms of the Peccei-Quinn scale f a and the initial misalignment angle θ i .Accounting for the remaining matter-dominated period between T a osc and T ϕ D so that and assuming no entropy production after T ϕ D , we arrive at the relic density estimate for misalignment axions: (2.31)

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 [60].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 ≡ m χ /T χ , can be approximated as [35,58] 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 for both the modulus and axion (i.e.α ∈ {ϕ, 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 [35,58] where we again use x ≡ m χ /T χ .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 [57].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 χ through a comparison of ρ χ to n χ m χ 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 temperature 6 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 The initial energy density of the modulus is given by where we take ϕ 0 = m P , as expected on dimensional grounds [42,43] and was found in an explicit calculation for Kähler inflation [44].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 ⟨σv⟩ χ ρ χ (T R )/m χ 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 ρ where we define c through the expected decay rate formula [18,56,59] (2.38) The number of relativistic degrees of freedom, g * (T ), is fitted from data used in [35,59].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. [37]  ) the modulus has an effective coupling parameter c that is compatible with explicit models. 7We 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.

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 In the absence of anisotropic stress for the fluid sources, we have Φ = Ψ.Working in momentum space and suppressing the wave number, the time-time component of the perturbed Einstein 7 We take c = 27.5 throughout this work based on results presented in [18].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 Γ ϕ is proportional to c, we will see that raising c significantly leads to a larger T ϕ D and thus makes our results even more constrained.
where the prime denotes derivatives with respect to number of e-folds N , and v (α) , δρ (α) , δp (α) are scalar velocity, density and pressure perturbations for each fluid, respectively.Introducing fractional density perturbations δ (α) ≡ δρ (α) /ρ (α) and defining the velocity perturbation for each fluid as θ (α) = a −1 ∇ 2 v (α) , the continuity equations in momentum space are given by Similarly, the equations for velocity perturbations are 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 [31], we set the modes after modulus domination when all modes of interest are super-Hubble, k/aH → 0. As discussed in [31], isocurvature constraints do not apply to the model we are considering here and so we are interested in strictly adiabatic initial conditions for the multifluid perturbations 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 ∼ H −1 ≫ m −1 ϕ .An approximate solution is then (3.12) where we choose initial values so that ρ ϕ 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 Γ ϕ /H ≪ 1, it follows that This relation differs from the standard relations due to the presence of decays and entropy production.Taking the super-horizon limit k/aH → 0 of Eq. (3.2) in a modulus dominated universe, we have where we used conservation of Φ on super-Hubble scales.Since ρ ϕ ≃ 3H 2 m 2 p during modulus domination, Eq. (3.16) implies the following initial condition for long wavelength gravitational perturbations δ (0) ϕ = −2Φ 0 and it follows from Eq. (3.11) that δ Finally, because scalar velocity perturbations quickly decay outside of the Hubble radius their initial values can be ignored.

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 −1 rh , 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 Φ = Φ 0 for both super and sub-Hubble scales.Using this, we can rewrite Eq. (3.7) as 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 which confirms that long wavelength vector modes are unimportant.From Eq. (3.18), we can derive the evolution of the modulus perturbation δ ϕ during the modulus dominated era subject to the initial condition δ (0) ϕ = −2Φ 0 .Noting that until the time of reheating we have Γ ϕ /H ≪ 1 and Φ is constant, we can rewrite Eq. (3.3) as and using the result in Eq. (3.18), we integrate to find 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 [31] where by using the background approximation Eqs.(3.12) -(3.14), we have time-dependent coefficients scaling as 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, It was shown in [31] 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 [31] Similarly, using Eq.(3.20) and Eq.(3.25) in Eq. (3.5) and remembering that the background coefficients are constants, one finds which again is valid on both super-Hubble and sub-Hubble scales, and we have used the initial conditions δ Given the solutions for the modulus and WIMP perturbations, we can solve for the radiation fluid perturbations.The equation for radiation perturbations is [31] δ ′′ r + 2A − where A ≡ A 3 + A 4 and α = 2 3 and the source term is given by A solution to Eq. (3.27) is which can be further simplified by integration by parts.We thus obtain In the absence of decays and annihilation (corresponding to A = α = 0 above), the exact solution to Eq. (3.27) simplifies significantly for all scales: 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 [31] for a more detailed analysis.

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 As the modulus decays become important, the constant scaling from the modulus dominated phase is no longer valid and we will perform analysis numerically following the approach in [31].
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 time 8 .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) [31,61], 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, 9 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 [31].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 [62].In the absence of WIMP annihilation (i.e. in the branching scenario), this result provides the seeds for further growth of structure [61], whereas when annihilation is important they act to further reduce the strength of the perturbations following reheating [31].
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 [31,61].
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 ∼ 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 freestreaming scale is determined by an integral involving the average velocity following decoupling until today t 0 as 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 where When the modulus decays to WIMPs, the WIMPs will thermalize if their scattering rate exceeds the Hubble rate.In this case, This leads to the washout 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 −1 f s .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 and the WIMPs will decouple prior to reheating.
There will typically be a continuous momentum spectrum and the average momentum at reheating ⟨p R ⟩ is model dependent.However, we can gain insight by considering when the DM resulting from decay is primarily non-relativistic ⟨p R ⟩ ≪ m χ .In this case one finds [31,61] whereas for the relativistic case ⟨p R ⟩ ≫ m χ the result is where the subscript 'eq' denotes the time of radiation and matter equality, and 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

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): where H(T a osc ) is the Hubble scale when axion oscillations begin.We have also defined our scale factor normalization to be a(T a osc ) ≡ 1.In [37], 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 δ a ∼ 1 by the end of the EMD phase.In this work, we also look at cases where δ a = 10 −2 and δ 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 [37] for a mode with wavelength k −1 entering the horizon at some temperature T a osc ≳ T ≳ T ϕ D .For modes which are subhorizon when axion oscillations begin, i.e. for k > H(T a osc ), 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 a osc ), sufficient growth of the perturbations may be difficult.We therefore focus on the modes which enter the horizon close to T a osc so that k ∼ H(T a osc ), 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: where −Φ 0 ∼ 10 −4 and we have the conditions a .Thus, for a given T ϕ D 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.We show in Fig. 3 this bound on the PQ scale for each δ a ∈ {1, 10 −2 , 10 −3 }.Here, we note that imposing larger values of δ a for a given f a also forces T ϕ D to be closer to the BBN limit.For example, for f a = 10 9 GeV, T ϕ D is constrained to be below 33 MeV for δ a = 1 while it may be as large as 240 MeV for δ a = 10 −2 and 780 MeV for δ a = 10 −3 .These bounds correspond to m ϕ ≲ 2.4 × 10 5 GeV, m ϕ ≲ 1.2 × 10 6 GeV, and m ϕ ≲ 2.8 × 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 δ a require longer periods of matter domination, and therefore lower T ϕ D .Conversely, for a given δ a , increasing f a (and thus lowering the scale of axion oscillations) requires lowering T ϕ D 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 Both qualitative behaviors described above can be understood from the parametric scaling in Eq. (4.4).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 a osc , Eq. (4.2) can be rewritten as where A 2 0 is given by Eq. (2.29).We show this estimate in Fig. 4. Here, we adopt T ϕ D and corresponding T a osc values from the bounds on (T ϕ D , f a ) pairs for the various δ a which are shown D allowed by smaller δ a require an increase the axion oscillation temperature, as can be seen from Eq. (2.27).This causes a reduction in m a (T a osc ) which, as can be seen in the figure, can compete against the previously mentioned reduction in m a (T ϕ D ).We also note that for δ a = 1, we do not display f a ≳ 10 16 GeV as the minicluster formation criteria can only be met for T ϕ D 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 δ a , it is possible that much larger superhorizon modes k < H(T a osc ) 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.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 for an untuned θ i = 2. Purple regions are consistent with the minicluster formation criterion δ a = 1, dark blue with δ a = 10 −2 , and light blue with δ a = 10 −3 .Orange regions correspond to even lower values of δ a but still have some overlap with an EMDE.Red regions, however, have T a osc < T ϕ D 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 > T a osc , 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 ≳ O(10 13 ) GeV.For lower values of f a , the horizontal dashed line displaying Ω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 ≲ O(100) MeV.As expected, the regions compatible with minicluster formation become smaller the larger δ a is.For δ a = 1, for example, we are forced into the purple regions which only intersect the relic density line when ϕ decays very close to BBN and f a ∼ O(10 16 ) GeV.Additionally, if T ϕ D is further required to be above ∼ 4 MeV for compatibility with BBN as suggested in [54,55], our results become even more constrained so that only f a ≲ O(10 14 ) GeV is allowed and the total relic density cannot be satisfied with the choice θ 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 means shorter periods of matter domination and hence smaller values of δ a ; thus, purple regions transition to dark blue, then light blue, and finally orange.On the other hand, for a given T ϕ D increasing f a also lowers δ 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, θ i → π, the axion relic abundance can be enhanced by a potentially large amount due to anharmonic effects [51].We display the axion relic density for a tuned θ 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 θ 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 ≲ O(200) MeV and f a ∼ 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.

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 δ a ∼ O(10 −3 − 1) for values of f a ≲ 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 δ 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 δ a ∼ 1, one ultimately runs out of parameter space even for large f a , as the BBN limit on T ϕ D 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: 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 ≲ m χ and the equilibrium density ρ χ becomes Boltzmann-suppressed, the precise value of ⟨σv⟩ χ 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 ≳ s 0 .Here we focus on the case of a purely s-wave annihilation cross section so that ⟨σv⟩ χ = 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 [63], MAGIC [64], and Fermi-LAT [65,66], which have constrained the T → 0 limit of ⟨σv⟩ χ for several different production channels.Here we take a model-independent approach and adopt the limits from the τ + τ − channel shown in [65].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.
Firstly, the Fermi-LAT upper bounds on the annihilation cross section to τ + τ − 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 ⟨σv⟩ multi χ in a multi-component scenario can be related to the constraint ⟨σv⟩ single χ in a WIMP-only scenario by [67] where ξ ≤ 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).
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 ϕ D ∼ O(1 − 10) MeV.We first understand the implications of 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).
this for WIMPs in the annihilation scenario, where the relic density is given by Eq. (2.24).For low T ϕ D and T f ∼ m χ /20 ∼ 5 GeV (for a 100 GeV WIMP), we have 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 ϕ→χ = 5% so that we are indeed in the annihilation scenario.We also set m ϕ = 4.1 × 10 5 GeV so that T ϕ D ≃ 51 MeV, roughly corresponding to the upper limit allowed by δ 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 ξ) mitigates the tension somewhat, but is generally unable to overcome the effect of the large value of T f /T ϕ 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., as in Fig. 7.We have also fixed B ϕ→χ = 5%.
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 (5.6) Thus, for T ϕ D ∼ O(1 − 50) MeV required for axion minicluster formation the branching of the modulus to WIMPs falls in the range B ϕ→χ ∼ 10 −4 − 10 −7 depending on the precise magnitude of m χ .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 ϕ→χ ∼ 10 −4 − 10 −6 , the annihilation scenario is also operational beginning at T ϕ D ≳ 0.2 GeV for small WIMP masses.This results from m χ being sufficiently low so that the late WIMP freeze-out temperature (close to T ϕ D ) results in ρ χ being large enough to produce WIMPs through thermal interactions.For larger m χ , the WIMP freeze-out temperature is then larger than T ϕ D 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 ϕ D that is compatible with axion minicluster formation, as obtained from Fig. 5.We see that for the δ a = 1 bound, the total DM relic density can indeed be saturated for the displayed branching ratios B ϕ→χ ∼ 10 −4 − 10 −6 for most of the range of m χ probed by Fermi-LAT.However as δ a decreases the parameter space becomes more restricted, requiring heavier WIMPs with m χ ≳ 1 − 2 TeV to have B ϕ→χ ≲ 10 −7 .Many currently existing models involving moduli and WIMPs predict B ϕ→χ ∼ 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 ϕ→χ .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.

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 osc > T ϕ D .Secondly, the axion minicluster formation condition Eq. ( 4.3) requires H(T a osc ) ≫ H(T ϕ D ).However, since H(T ϕ D ) ∼ Γ ϕ ∝ (T ϕ D ) 2 and H(T a osc ) ∼ m a (T a osc ) ∝ f −1 a , these conditions effectively force T ϕ 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 ϕ→χ ≲ 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 ϕ→χ for minicluster formation may also predict WIMP DM substructuresespecially 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 → 0 limit on ⟨σv⟩ χ 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.

Figure 3 :
Figure 3: The maximal value of f a which can be expected to form miniclusters for a given T ϕ D .The condition for minicluster formation is taken to be δ a ∈ {1, 10 −2 , 10 −3 } for the solid, dashed, and dot-dashed curves respectively.

Figure 4 :
Figure 4: The estimated minicluster mass for modes with k ∼ H(T a osc ) at the time of formation, T ϕ D , for a given f a .Superhorizion modes which have k ≳ H(T a osc ) may produce larger miniclusters than those displayed.The condition for minicluster formation is taken to be δ a ∈ {1, 10 −2 , 10 −3 } for the solid, dashed, and dot-dashed curves respectively.

Figure 5 :
Figure 5: Axion relic density versus T ϕ D where θ i = 2.For our benchmark, these values of T ϕ D correspond to m ϕ ∈ [2.4 × 10 4 , 5 × 10 6 ] GeV.Purple regions are consistent with δ a = 1, dark blue with δ a = 10 −2 , and light blue with δ a = 10 −3 .Orange regions correspond to even lower values of δ a but still have some overlap with an EMDE.Red regions have no overlap with an EMDE.

Figure 6 :
Figure 6: Axion relic density versus T ϕ D where θ i = 3.113.For our benchmark, these values of T ϕ D correspond to m ϕ ∈ [2.4 × 10 4 , 5 × 10 6 ] GeV.Purple regions are consistent with δ a = 1, dark blue with δ a = 10 −2 , and light blue with δ a = 10 −3 .Orange regions correspond to even lower values of δ a but still have some overlap with an EMDE.Red regions have no overlap with an EMDE.

Figure 9 :
Figure 9: Branching Scenario: We display contours of the branching of the modulus to WIMPs B ϕ→χ for which the observed DM relic density is obtained, assuming Fermi-LAT constraints with ξ = 1.Here, the solid boundaries of contours show where the overall DM relic density is saturated exclusively by WIMPs, while shaded regions have Ω χ h 2 < 0.12.For T ϕ D ≳ 0.2 GeV, the scenario transitions to the annihilation scenario for the values of B ϕ→χ shown.The black vertical lines indicate the highest T ϕ D for which axion minicluster formation is expected to occur for δ a ∈ {1, 10 −2 , 10 −3 }.