Multiverse Dark Matter: SUSY or Axions

The observed values of the cosmological constant {\it and} the abundance of Dark Matter (DM) can be successfully understood, using certain measures, by imposing the anthropic requirement that density perturbations go non-linear and virialize to form halos. This requires a probability distribution favoring low amounts of DM, i.e. low values of the PQ scale $f$ for the QCD axion and low values of the superpartner mass scale $\tilde{m}$ for LSP thermal relics. In theories with independent scanning of multiple DM components, there is a high probability for DM to be dominated by a single component. For example, with independent scanning of $f$ and $\tilde{m}$, TeV-scale LSP DM and an axion solution to the strong CP problem are unlikely to coexist. With thermal LSP DM, the scheme allows an understanding of a Little SUSY Hierarchy with multi-TeV superpartners. Alternatively, with axion DM, PQ breaking before (after) inflation leads to $f$ typically below (below) the projected range of the current ADMX experiment of $f = (3 - 30) \times 10^{11}$ GeV, providing strong motivation to develop experimental techniques for probing lower $f$.


Introduction and Overview
Our universe evolved from an era of radiation domination to one of matter domination and is currently transitioning to one of cosmological constant domination. Thus, over at least ten decades of temperature, T , the energy density of the universe involves two unknown parameters ρ = Λ + ρ m + ρ rad ∼ Λ + ξT 3 + T 4 (1) where, in the final form, numerical factors are omitted. The matter abundance is parameterized by ξ, the ratio of the matter energy density to the photon number density ξ = ρ m /n γ , and has both baryon and DM components ξ = ξ D + ξ B . The simplicity of this result should not mask the presence of two key theoretical problems

• The Cosmological Constant
Why is Λ so small?
The cosmological constant problem has a simple anthropic solution [1]: typical cosmological density perturbations do not become non-linear unless where Q ∼ 2 × 10 −5 is the magnitude of the primordial perturbations as they enter the horizon and f = 1 for ξ D ξ B . The absence of virialized matter halos leads to a dilute gas of inflating particles without any large scale structure or observers. While it is convincing that (2) is a necessary requirement for observers, the expected multiverse probability distribution for Λ, dP ∝ dΛ, leads to Λ typically two to three orders of magnitude larger than observed [2]. Furthermore, this condition does not directly address the Why Now problem; there could be an arbitrary long delay between the formation of halos and the appearance of observers.
The multiverse produced by eternal inflation is infinite; hence, when computing probabilities it is necessary to regulate divergences. In the Causal Patch measure [3], which removes the divergence by looking at finite regions around geodesics, the number of observers is diluted by inflation unless the time of Λ-domination occurs after the era at which observers occur, t Λ > t obs [4,5,6]. Since the dilution is exponential, this measure effectively forces the anthropic constraint where G is Newton's constant. Remarkably, this solves both the Cosmological Constant and the Why Now problems [4,5,6]: the most probable observed universes are those close to the boundary (3), having t Λ ∼ t obs . Today we are observing the onset of Λ-domination, so that the prediction of Λ from (3) is more successful than from (2). Hence, in this paper we assume a measure that leads to (3); currently the only known such measure is the Causal Patch measure.
A prediction for the cosmological constant, whether by (2) or (3), follows if the only scanning parameter in the multiverse is Λ, and is apparently destroyed if other relevant parameters, (Q, ξ D , ξ B , G, t obs ), scan. To prevent runaways there must be additional environmental selection. For example, requirements of galactic halo formation and cooling were included in multiparameter scans in [7,8]. In a more general scanning of parameters, one hopes to discover that our universe lies at the tip of a multi-dimensional cone formed by anthropic boundaries from a variety of cosmological and particle physics requirements [9]. For example, the virialization constraint (2) is a surface in this multi-dimensional space; a slice through the parameter space at fixed (Q, ξ D , ξ B ) leads to the original interpretation [1] of (2) as a bound on Λ.
In general, the multiverse probability distribution depends on n scanning parameters, and can be integrated over m parameters to yield an effective distribution for the remaining n − m parameters. Differing subspaces are of interest for different problems. For the evolution of our universe the subspace (ξ D , ξ B , Λ) is key, since these are the parameters that enter (1). For simplicity, in this paper we set the baryon density to its observed value, ξ B = ξ B0 , and scan in the 2d subspace (ξ D , Λ).
We take a probability distribution where p(ξ D ) Λ is the a priori distribution, suitably marginalized over other scanning parameters, while n obs is the environmental weighting factor that we take to contain three contributions n obs (ξ D , Λ) = n Λ meas (Λ) n ξ D meas (ξ D ) n vir (Λ, ξ D ). (5) n Λ meas contains the exponential dilution factor from the measure that imposes (3), while results from the dilution of baryonic observers for ξ D > ξ B [10,11] implied by such measures, in particular by the causal patch measure. Finally, n vir (Λ, ξ D ) is the environmental weighting factor that follows from requiring density perturbations to go non-linear and virialize. In much of the (ξ D , Λ) plane the probability is exponentially suppressed; at low ξ D from the tail of the Gaussian distribution of the initial density perturbations, and at large Λ from dilution of observers by inflation. The onset of this exponential behavior can be approximated by catastrophic boundaries, and in Figure 1 we show the virialization and observer-dilution boundaries from the exponential behavior of n vir (Λ, ξ D ) and n Λ meas (Λ). These boundaries correspond to (2) and (3). The transition at these boundaries is quite sudden; as the shaded regions of Figure 1 are entered, the probability of observers is exponentially reduced.
On the other hand we assume that the a priori distribution p(ξ D )Λ varies much more slowly and, in the region of interest, we take p(ξ D ) to be a power law. In the limit of taking sharp catastrophic boundaries, the probability distribution in the unshaded observer region is f (ξ D , Λ) = p(ξ D )Λ n ξ D meas (ξ D ). If the probability force (∂ ln f /∂ ln ξ D , ∂ ln f /∂ ln Λ) points towards the catastrophic boundaries, as shown by the red arrow in Figure 1, typical observers should lie near the tip of the cone formed by the intersection of the two boundaries, as is the case for our universe shown by the red star. Shaded regions with solid line borders have exponential suppression of observers by at least 1%, while the dashed lines are for 5% suppression, as discussed in more detail in section 2. A multiverse probability force, illustrated by the red arrow, allows our universe, the red star, to be typical. The component of the force to large Λ is predicted, while that to low ξ D is assumed.
At fixed ξ D = ξ D0 , our universe is a factor of 10 2 (10 3 ) in Λ from the most probable universes with halos of mass M = 10 12 (10 7 ) M [2,4]. Similarly, from Figure 1 our universe is a factor of 10 2 (10 3 ) in Λ from the virialization boundary for halos of mass M = 10 12 (10 7 ) M . However, as the virialization boundary is so steep, reflecting the fourth power of ξ in (2), at fixed Λ = Λ 0 , our universe is only a factor 3-5 from the boundary in ξ D . Hence, it is reasonable to take the view that Λ is determined by observer dilution and ξ D by virialization. In this 2d plane we simultaneously make use of (2) and (3) to understand both the size of the cosmological constant and the amount of DM We stress that the virialization boundary, used in [1] to predict a non-zero cosmological constant, is here used to determine the abundance of DM.
We note that ξ D is the relevant variable to describe the abundance of DM; it is the quantity computed directly in any theory of DM genesis. For example, models of WIMP DM characterized by a single mass scale m lead to a freeze-out abundance ξ D ∝ m 2 . In this example, the DM abundance in our universe corresponds to a value for m approximately a factor of 2 from the virialization boundary.
Other boundaries, in particular the requirement of halo cooling [7,8], could play an important role in the (ξ D , Λ) plane. However, the physics of this boundary is more complicated than for the virialization and observer-dilution boundaries, involving certain details of atomic and molecular physics, merger statistics and shocks, and has power dependence on the halo mass. Hence we neglect halo cooling in this paper; inclusion would not significantly change our results.
In section 2 we review the virialization boundary. We stress that starting with our universe and decreasing ξ D leads rapidly to the loss of large scale structure. In our view, virialization forces our universe to contain DM. Of course, this view requires a distribution p(ξ D ) favoring low ξ D -there is a cost in the multiverse to produce DM. In the rest of the paper we explore the consequences of this important requirement. 1 In section 3 an effective distribution for ξ D is obtained by integrating over the cosmological constant. For conventional LSP thermal freeze-out, the effective distribution for the superpartner mass scalem should favor low values and, depending on the model, could lead to the TeV scale. For QCD axion DM the effective distribution should favor a low PQ breaking scale, frequently leading to f ∼ (10 10 − 10 11 ) GeV.
In section 4 we show that a probability distribution favoring low ξ D leads to a strong expectation that DM is dominated by a single component. 2 To be close to the virialization boundary we argue that all components of DM must have distributions favoring low values, and the dominant component will be the one with weakest distribution, with the others highly subdominant. For a supersymmetric theory with both a thermal freeze-out LSP and an axion this would mean eitherm ∼ TeV and f (10 10 − 10 11 ) GeV orm TeV and f ∼ (10 10 − 10 11 ) GeV; both are experimentally excluded.
LSP DM in supersymmetric theories withm scanning is considered in more detail in section 5. In the case that weak-scale superpartners lead to too little DM, as would happen with a wino or Higgsino LSP, the scale of superpartners will necessarily be lifted to allow sufficient DM for virialization to occur, thereby generating a Little Susy Hierarchy with multi-TeV superpartners.
In section 6 we study QCD axion DM in more detail. We begin by elucidating the parameter space of this cosmology, including both Pre-and Post-Inflation cases. We determine which power law probability distributions for f lead to the observed proximity of our universe to the virialization boundary, and obtain statistical predictions for the value of f in our universe in the Pre-Inflation cosmology. We find f is frequently (always) below the expected reach of the ADMX experiment for Pre-(Post-) Inflation cosmologies. New experimental techniques to probe this low f window are urgently needed. For our universe to be typical with a distribution favoring low f , thermal axion universes with f ∼ 10 5±1 GeV must be catastrophic for observers. We study thermal axion universes and conclude that observers are likely to be suppressed.

Halo Virialization: the requirement of Dark Matter
Halo Virialization on comoving scale λ occurred when the matter density perturbation δ m (λ) went non-linear. Taking a wide view of the history of our universe this is a recent phenomenon. Perturbations on the scale of our galaxy typically went non-linear at a redshift of a few. As we look at other regions of the multiverse, we find that virialization is far from guaranteed. The vast majority of the multiverse has values of the cosmological constant, Λ, far too large to allow relevant perturbations to go non-linear. We begin by considering universes where the amount of DM varies but all other fundamental parameters, including Λ, are fixed at their observed values. Later in this section we also scan Λ and consider the virialization boundary in the (ξ D , Λ) plane.
ξ B,D are defined to be the ratio of the baryon and DM energy densities to the number density of photons and are determined by baryogenesis and DM genesis. Varying ξ D induces a variation in T e,Λ , the temperatures of matter-radiation and matter-Λ equality, as well as the time-temperature relation. 3 Ignoring all numerical factors, the parametric form for the energy density is ρ ∼ T 4 + ξT 3 + Λ, with terms for radiation, matter and the cosmological constant, and ξ = ξ B + ξ D .
Hence T e ∼ ξ and T Λ ∼ (Λ/ξ) 1/3 . The evolution of density perturbations is well-known and has been studied at a high level of accuracy. In Figure 2, using a fitting formula of [16], we plot the ratio δ m /δ m0 , at an era that is well after recombination and matter-radiation equality but still in the perturbative linear regime, against ξ D /ξ D0 . (The subscript "0" refers to values in our universe.) The four curves are for comoving scales of k = (0.1 − 1)Mpc −1 , corresponding to baryonic masses of M ∼ (10 15 − 10 12 )M , where M is the mass of our sun. The most striking feature is that a reduction in the amount of DM causes a very substantial drop in the density perturbation and, for masses of our galaxy and smaller, this reduction becomes extreme at very low DM abundances.
This behavior can be understood from the rough analytical approximation where f B,D = ξ B,D /(ξ D + ξ B ) are the fractions of matter contained in baryon and dark species, and the density perturbations are taken to enter the horizon with size Q 0 = 2×10 −5 , independent (which assumes the linear regime). The two terms in (9), proportional to f B and f D , originate from baryon and DM perturbations and evolve very differently. On mass scales less than M S the baryon perturbations are reduced by photon diffusion (Silk) damping near the era of recombination, M S0 ∼ 10 16 M . On the other hand the DM is not coupled to the photons and is able to grow logarithmically during the radiation dominated era by a factor G rad . After recombination, the baryons are released from the coupling to photons and the perturbations of baryons and DM can be combined into a single matter perturbation, δ m of Eq. (9). In our universe the baryon perturbations on scales relevant for virialization are strongly Silk damped, but this is almost irrelevant as after recombination the baryons fall into the potential of the DM perturbations, which are large since f D G rad f B . Returning to Figure 2, the lower curve is for halo mass M ∼ 10 12 M , sufficiently low that the baryon contribution to the perturbation is negligible and the resulting perturbation drops linearly with f D . The curves at successively larger M eventually reach a plateau where the baryon contribution, even though damped, dominates the dark perturbation at low f D .
As ξ D is reduced below ξ D0 the reduction in δ m is immediate. Even though f D at first stays close to unity, the matter dominated growth factor T e /T Λ scales as ξ 4/3 D . Once the DM density drops below that of baryons, T e /T Λ becomes constant, but f D scales as ξ D and drops below unity. The second term of (9) is then reduced because, when the baryons are released from the photon coupling after recombination, there is only a small potential from the DM for them to fall into. The numerical solution shown in Figure 2 is smooth at ξ D = ξ B0 , but asymptotically the slope varies from 4/3 to 1, for perturbations where the baryon component is Silk damped.
From Figure 2 we see that when ξ D /ξ D0 ∼ 1/5 the matter perturbation on galactic scales and smaller is already suppressed by an order of magnitude. This means that typical perturbations on these scales never go non-linear as their growth is halted by the onset of Λ domination. On rare occasions the perturbations may fluctuate to larger values, allowing regions of growth, but this becomes exponentially rare as ξ D drops further.
We obtain a virialization boundary, shown in Figure 1, by fitting the curves in Figure 2 of [4] by an exponential. With good accuracy the n vir (Λ, ξ D ) factor depends only on the halo mass and the value of the matter density perturbation δ m (t Λ ) at the time of cosmological constant domination. We consider two values for the halo mass M = 10 (7,12) M and find The solid (dashed) lines in Figure 1 correspond to those points in the (ξ D , Λ) plane where the exponential factor in (11) drops to 1% (5%). Thus the excluded regions contain 1% (5%) of observers. For example, for Λ = Λ 0 and M = 10 (7,12) M we find the 1% boundary to be ξ 1% = (0.10, 0.23)ξ 0 and the 5% boundary to be ξ 5% = (0.12, 0.27)ξ 0 . Figure 1 also shows the observer dilution boundary. We define this by extracting the exponential behavior of n Λ meas (Λ) from [5]: Λ obs is related to the time t obs at which observers occur by The solid (dashed) horizontal line in Figure 1 corresponds to a value of the exponential factor in (12) of 1% (5%). Baryonic perturbations on scales larger than M S are not subject to Silk damping, and therefore could be relevant for large scale structure even if ξ D = 0. If such perturbations go non-linear it is possible that a top-down scheme for large scale structure could lead to observers. The difficulty is estimating the probability of such observers relative to ones in universes similar to our own. Since we assume a probability distribution for ξ D that decreases as ξ D drops below ξ D0 , there will be a gain in probability from reducing ξ D well below ξ B0 . However this can be offset by the reduction in Λ necessary to allow such large baryonic perturbations to go non-linear. With ξ D irrelevantly small, T e drops by about a factor 6 from T e0 , and such large baryonic perturbations don't grow significantly until the matter dominated era is reached. We estimate that such perturbations go non-linear only if T Λ is reduced by ∼ 10 2 relative to our universe and the reduction in Λ to achieve this costs ∼ 10 −6 in probability. Halo and star formation is entirely different in a top-down scheme, so that the corresponding observer weighting factors will be very different and are unknown.

An Effective Distribution for ξ D
For the rest of the paper it is convenient to perform an approximate integral of the probability distribution over Λ so that we can focus our discussion on the remaining effective probability distribution for ξ D . We define the observer region to be where the exponentials in n vir (10 12 M ) and n Λ meas , shown in (11) and (12), are larger than ∼ 50%. Thus the observer region is similar but slightly smaller than the unshaded region in Figure 1, where the boundaries were defined by these exponentials reaching 1%. We define the observer dilution and virialization boundaries of the observer region to be so that they intersect at (ξ c , Λ c ) ∼ (0.5 ξ D0 , Λ 0 ). Approximating n vir and n Λ meas as θ functions at the boundaries (14), in the observer region the probability distribution of (4, 5, 6) becomes Here we use the value of ξ B observed in our universe as we are not scanning over the baryon density.
Integrating over Λ gives an effective distribution for the single parameter ξ D where the last factor arises because the maximum value of Λ in the observer region depends on ξ D . We assume that the (ξ D /ξ c ) 4 factor makes the region ξ D < ξ c sufficiently improbable that there is little error in taking the 1d observer region to be This assumption would only fail if the distribution p(ξ D ) favored low ξ D so strongly that it overcompensates the (ξ D /ξ c ) 4 factor. However, this would cause runaway along the virialization boundary in Figure 1 to both low ξ D and low Λ, destroying the successful understanding of the Why Now problem. For less extreme p(ξ D ) distributions, the Why Now problem is solved, since the typical value of Λ is Λ c ∼ Λ 0 . In the rest of the paper we fix Λ = Λ 0 and perform a 1d scan over ξ D with distribution (16) subject to (17).

Single Component Dark Matter
The arguments of previous sections are independent of the nature of DM. In this section we show that if the abundance of DM is explained by the virialization boundary it is dominantly composed of a single component. Suppose that there are many components, ξ D = Σ i ξ i , with sufficiently many scanning parameters that the energy densities, ξ i , scan independently with prior dP = Π i P i (ξ i ) d ln ξ i . In our universe, if some species has ξ i ξ c density perturbations are easily able to go non-linear, and we are far from the virialization boundary. Hence, for the abundance of DM to be explained by proximity to the virialization boundary, all ξ i should be of order ξ c or smaller. Is it possible that more than one component has ξ i ∼ ξ c , with distributions favoring low ξ i ? In general this possibility is extremely unlikely. The component with the weakest distribution towards the boundary will have an abundance close to ξ c , while the other components will typically have much smaller abundances, as illustrated in Figure 3. Hence, if the virialization boundary is the correct explanation for the DM density, DM is strongly dominated by a single component. 4 How small might we expect the sub-dominant components to be? In generic theories of DM, ξ i can vary over many orders of magnitude. A sub-dominant species must have the gradient dP i /dξ i negative at ξ i = ξ c so that ξ i will runaway to low some low value ξ min i where the sign of dP i /dξ i changes. It would be accidental for this to happen anywhere near ξ c , which is determined by the physics of virialization and is independent of the prior distribution. Hence, we expect ξ min i to be less than ξ c by at least a few orders of magnitude, as illustrated by the stars in Figure  3.
How improbable is multi-component DM? Consider a two-component model with densities ξ 1,2 scanning independently with an effective distribution with n 1,2 < 0 and C a normalization constant. We take (18) to be valid over a wide range of ξ 1,2 , including the region of the virialization boundary near ξ 1 ∼ ξ 2 , and down to near ξ min 1,2 , where it breaks down and the gradient of the distribution changes sign. Let P multi be the probability for ξ 1,2 ∼ ξ c integrated over a region with ∆ξ 1,2 ∼ ξ 1,2 . Similarly let P single 1 be the probability for (ξ 1 , ξ 2 ) ∼ (ξ c , ξ min 2 ) integrated over a region with ∆ξ 1,2 ∼ ξ 1,2 . We find In going from single component DM with (ξ 1 , ξ 2 ) ∼ (ξ c , ξ min 2 ) to multi-component DM with (ξ 1 , ξ 2 ) ∼ (ξ c , ξ c ) one loses in probability from the ξ n 2 2 factor of (18) by changing ξ 2 from its minimal value to ξ c . The result for P multi /P single 2 is obtained by interchanging 1 ↔ 2. For The blue arrows represent one choice for ∇P i which leads to runaway along the virialization boundary to small ξ 1 , as shown by the blue star, and the red arrows represent another choice that leads to runaway to small ξ 2 and the red star.
example, for |n i | = N and ξ min i = 10 −3 ξ c , multi-component DM is less probable than single component by a factor of 10 3N . Consider applying these general arguments to supersymmetric theories that contain a QCD axion, and assume a cosmology that leads to both misalignment axion DM, ξ a , and thermal LSP freezeout DM, ξ LSP . Further assume that the multiverse distribution functions for the axion decay constant, f , and the scale of supersymmetry breaking,m, favor low values, so that the size of ξ D = ξ a + ξ LSP will be explained by proximity to the virialization boundary (since ξ a ∝ f and ξ LSP ∝m 2 ). The relative strengths of the distributions for f andm are unknown, so we do not know which is the dominant component with density near ξ c and which the sub-dominant component with density at least a few orders of magnitude below ξ c . If axions are dominant f ∼ 10 11 GeV, while if LSPs are dominantm ∼ TeV, (these are rough order of magnitude estimates, which have theoretical uncertainties and model dependencies, but are sufficient for our purposes). This then implies that if axions are dominantm is at least a few orders of magnitude below the TeV scale, while if LSPs are dominant f is at least a few orders of magnitude below 10 11 GeV. Both cases are observationally excluded -they do not describe our universe. Hence, in our scheme, with the abundance of DM explained by closeness to the virialization boundary, conventional TeV-scale LSP freezeout DM is inconsistent with the axion solution to the strong CP problem. To avoid this conclusion, ξ min a,LSP must both be accidentally close to ξ c , implying that the observed DM abundance is close to the peak of the 2d probability distribution.
Of course, theories can contain either stable TeV-scale LSPs or axions, and these both remain interesting DM candidates when the abundance is explained by the virialization boundary. We consider the dominant DM component to be the lightest superpartner in section 5 and the axion in section 6.

LSP Dark Matter and a Little SUSY Hierarchy
In the Standard Model the weak scale, v, is fine-tuned. Without a multiverse this suggests that some new physics, such as supersymmetry, should occur at the weak scale. However, in a multiverse such new physics is not required, since a finely-tuned weak scale may result from anthropic requirements [13,14]. Whether new physics is likely at the weak scale depends on probability distributions.
Consider a supersymmetric extension of the Standard Model with the overall mass scale of the superpartners scanning, but with fixed superpartner mass ratios. For convenience we definem to be the mass of the Lightest SuperPartner (LSP), which we assume is cosmologically stable and has a thermal freeze-out abundance ξ D ∝m 2 . Using (16,17) and assuming a prior distributionm n , we obtain an effective multiverse probability distribution The first factor in parenthesis is the measure factor of (6), withm 0 the value of the LSP mass that leads to ξ D = ξ D0 . The last factor in parenthesis arises from fine-tuning, to satisfy the anthropic electroweak symmetry breaking requirement, and form v becomes v 2 /m 2 . It is an additional factor in the environmental weighting quantity n obs , omitted in (5) since it is special to LSP DM. For simplicity we have assumed that the superpartners relevant for this tuning have mass close tom. As discussed in section 3, the argument of the theta function that approximates the virialization boundary corresponds to the weighting factor n vir = 0.5, leading tom c 0.7m 0 .
A crucial question is the size ofm 0 , which then determinesm c . So far our argument is quite general, with little dependence on the details of the supersymmetric model. However,m 0 is model dependent, depending on the nature of the LSP and its interactions. As is well-known, freeze-out frequently involves a mass scale somewhat larger than the weak scale. For example, m 0 1(3) TeV for pure Higgsino (wino) DM. We assumem 0 is larger than the weak scale and, for illustration, we takem 0 = 2.2 TeV, which givesm c = 1.5 TeV.
In Figure 4 the dashed curves show the probability distribution (20) for various n, with the θ function representing the virialization boundary ignored. For n < 0 the distributions are peaked at lowm. For 0 < n < 2 they are peaked near v, so that superpartners are expected at the weak scale. For 2 < n < 4 they are peaked nearm 0 / √ 6 ∼ 900 GeV, as illustrated by the dashed purple line. It is in this range of n that the multiverse predicts comparable amounts of dark and baryonic matter since the peak in the distribution is determined by the 1/(1 + ξ D /ξ B ) factor [11]. However, as shown by the shading in Figure 4, these distributions are Of course the detailed superpartner spectrum depends on the model, and may be compressed or split to varying degrees. We stress that a wide range of prior distributions, even those that strongly favor large values of supersymmetry breaking, lead to proximity to the virialization boundary and to a Little Hierarchy. 5 For n < 0 most universes havem v and in this region ξ D may have a different dependence onm. For example, if LSP annihilation occurs via virtual weak gauge bosons the annihilation cross section scales asm 2 /v 4 so that ξ D ∝m −2 , giving sufficient DM to allow large scale structure formation atm v. However, in many supersymmetric theories, such as the MSSM, it is not possible (or requires a costly fine-tune) to obtain v m, so these universes are excluded from anthropic requirements on the size of the weak scale. Hence, in such supersymmetric theories, with n < 0, one again has typical observers withm just abovem c .

Axion Dark Matter
The QCD axion [18,19] provides a minimal extension of the Standard Model (SM) that provides both a solution of the strong CP problem [20] and cosmological cold DM [21,22,23]. Over a period of several decades the importance of a QCD axion has steadily grown due to the absence of any other plausible solution to the strong CP problem, together with a succession of null results searching for weak-scale DM from both galactic DM and particle accelerators. The axion is a pseudo-Goldstone boson produced at scale f by the spontaneous breaking of a Peccei-Quinn (PQ) symmetry that carries a QCD anomaly.
We fix our notation by introducing the effective lagrangian we will use to discuss axion physics. The PQ symmetry is broken dominantly by a canonically normalized scalar field with vacuum structure φ = V e ia(x)/V , where a(x) is the axion field. At energies below the PQ breaking scale and the heavy SM fermion masses, but above the QCD confinement scale, the most general axion interactions are described by the following effective lagrangian at order 1/f The sum is over the light SM fermion fields ψ i = u, d, e, ν. Notice that (21) defines the constant f through the normalization of the GG term; it is related to the vev V by V = N f , where N is an integer determined by the color anomaly of the PQ symmetry, and is known as the domain wall number. QCD instantons and the PQ color anomaly generate a small zero-temperature axion potential ≈ Λ 4 QCD (1 − cos a/f ), leading to N equivalent vacua, and an axion mass that can be computed reliably using chiral perturbation theory

The Misalignment Mechanism and the Virialization Boundary
Since all axion interactions are suppressed by the scale f , for large values of f the axion is decoupled from the thermal bath in the early universe. The behavior of the classical axion field during these epochs is characterized by two phases. For H m a , the axion field is stuck at a random position in field space a(x) = a i . As soon as H m a the Hubble friction becomes irrelevant and the axion field starts to behave as in flat space oscillating around its equilibrium position, a = 0. Once the anharmonic corrections to the potential are small the energy density stored in these oscillations dilutes as non-relativistic matter with the Hubble flow and contributes as a cold component of the DM. This mechanism, which generates a contribution to the cosmological cold DM density through the coherent oscillations of a classical scalar field, is called the "misalignment mechanism".
In the "Pre-Inflation" cosmology, where the PQ phase transition occurs before or during inflation, initial fluctuations of the axion field on scales corresponding to our Hubble horizon can be neglected, so that the initial misalignment angle θ = a i /f is a constant. The axion DM energy density from misalignment has been computed to be [24,25] ρ Pre mis (f, θ) ≈ (2.0 × 10 −3 eV) 4 f 10 12 GeV m θ 2 F (θ) (23) with m 1.2. This result applies for f 10 15÷16 GeV, so that the axion field starts to oscillate at a temperature above Λ QCD , and includes all values of f of interest to us. The function F (θ) corrects for the anharmonicities of the potential, and has normalization F (0) = 1. This result has theoretical uncertainties arising from a non-perturbative instanton calculation of the temperature dependent axion mass [26]. Hence we introduce a corresponding uncertainty in the value of f that accounts for DM in the Pre-Inflation cosmology, and for illustration adopt the range Since misalignment axion production occurs near the QCD scale, the DM abundance of (23) is quite robust to non-standard cosmologies at higher temperatures. However, the abundance could be affected by changing the cosmology below the QCD scale; for example, entropy production would require larger values of f to explain the observed DM abundance.
In the "Post-Inflation" cosmology, where the PQ phase transition occurs after inflation, (23) has to be averaged over the possible values of θ where the initial misalignment angle is assumed to be uniformly distributed between -π and π, yielding If this were the only source of axion DM, the required value of f would be centered on 1.5 × 10 11 GeV. However, a network of axion strings is formed at the PQ phase transition and this network evolves, emitting axions. At the QCD phase transition each string becomes attached to N domain walls. Throughout our discussion of Post-Inflation axion cosmology we assume a domain wall number of unity, N = 1, so that there are no stable domain walls [27]. It is an interesting question whether N = 1 universes are anthropically disfavored. With N = 1, at the QCD phase transition strings become the boundaries of domain walls, which collapse and disappear, again radiating axions. Numerical simulations have studied the amount of axion cold DM arising from string evolution [28,29,30] and N = 1 domain wall collapse [31], reducing the value of f needed to account for DM and introducing further uncertainties in f . For illustration we adopt the range f Post = (10 10 − 10 11 ) GeV. (27) The variables distinguishing between the two previous regimes, whether to average or not over the values of θ, are the dynamics of inflation and reheating. In order for a i to be homogeneous   (17) for Pre-and Post-Inflation axion cosmologies respectively, with light shading corresponding to the theoretical uncertainties of (24) and (27).
within a Hubble patch the PQ phase transition has to have occurred already while the universe is inflating and the PQ symmetry must not be restored during reheating. These two facts require both T I ≡ H I /2π, the Gibbons-Hawking temperature of inflating de Sitter space, and T max , the maximal temperature reached by the universe during reheating, to be smaller than f . If inflation will insure that a i is constant within our Hubble patch. This possibility allows the so called "anthropic axion window" 6 with larger values of f but a small angle θ. 7 In Figure 5 the shaded region in the (f, θ) plane is excluded by the virialization boundary (17). In gray we plot the exclusion for the Pre-Inflation case in which the initial misalignment angle is constant in our Hubble patch. The thickness of the boundary of the excluded region shows the uncertainty from non-perturbative QCD corresponding to (24). The red shading corresponds to the exclusion for the Post-Inflation regime in which the misalignment angle is averaged over our Hubble patch. The theoretical uncertainty in this case (shown by the thick 6 The anthropic axion window [32] does not require a multiverse, but follows from conventional inflation leading to large regions having different θ when the PQ phase transition occurs before inflation. At large f an anthropic requirement that ξ D not be too large selects for small θ. This anthropic requirement is much less well understood than the requirement of virialization used in this paper. 7 An additional source of small scale fluctuation of the axion field are the quantum fluctuation which are imprinted on it during inflation. They determine a minimal value of the effective θ angle [33,34] θ min = H I 2πf ≈ 1.6 × 10 −4 H I 10 12 GeV and in turn a minimal contribution to ρ a in the Pre-Inflation case. boundary of the excluded region) also includes a contribution from the axion density from radiation from axion strings and domain walls, and corresponds to the range of (27). Figure 5 does not yet capture the full parameter space of axion DM as it misses the effect of T I and T max , the two parameters distinguishing between the Pre-and Post-Inflation scenarios. Their role is shown in Figure 6. In each of the four panels a different value for the energy scale of inflation, E I , is chosen. Its relation to the inflationary Hubble parameter H I and to the Gibbons-Hawking temperature is through the Friedmann equation In each panel the physical region is the one in which the temperature at the end of inflation, T max , is smaller than E I . The Pre-and Post-Inflation scenarios are separated by the boundary on which f = max(T I , T max ). The vertical part of the boundary has f = T I , and moves to lower f as E I is reduced from one panel to the next. The remaining part of the boundary is the straight line of unit slope, f = T max . The red (gray) shaded regions are excluded by the virialization boundary (17) in the Post-(Pre-) Inflation cosmology. In each case the lighter shading corresponds to the theoretical uncertainty corresponding to (27,24). For Pre-Inflation the position of the virialization boundary depends on θ and is shown for θ = 3 in all panels and for θ = 1 by the lightest gray shading in the lower two panels. A further constraint in the Pre-Inflation scenario arises from the generation of axion isocurvature perturbations during inflation [34]; the regions shaded blue in Figure 6 are observationally excluded for θ = 3. These perturbations grow rapidly with E I , and at large E I there may also be anthropic constraints.
To explain ξ D from the virialization boundary we take an effective distribution that favors low f . Figure 6 then shows that the Post-(Pre-) Inflation cosmology results for larger (smaller) values of E I and T max . For E I < 10 11 GeV, only the Pre-Inflation case is possible, while for E I > 10 14 GeV the description of our universe by Pre-Inflation cosmology becomes less probable.

Scanning over f and θ
In the previous sub-section we have shown the regions of axion parameter space excluded by the anthropic requirement that DM density perturbations go non-linear, allowing virialization and halo formation. As discussed in the introduction, in this paper we explore the possibility that the observed DM energy density is determined by this virialization boundary, taking the cosmological constant to be anthropically constrained by observer dilution, as illustrated in Figure 1.
We now investigate allowing f to scan over the multiverse according to a given prior probability distribution dP ∝ p(f )d ln f.
We stress that this is rarely studied, and is not the case in the conventional "anthropic window" where only θ scans. In the Pre-Inflation case the initial value of the misalignment angle is an additional parameter, with a flat probability distribution between −π and π. From Figure 1 we  In the Post-Inflation regions the DM abundance is independent of θ, while in the Pre-Inflation region it is proportional to θ 2 F (θ). The dark red (dark gray) shading shows regions of Post-(Pre-) Inflation cosmology excluded by the virialization constraint, with the lighter shading corresponding to the uncertainty of eq. 27 (24). Pre-Inflation shading is shown for θ = 3, and in the lower two panels the light gray region with dashed boundaries shows the uncertainty band for θ = 1. Blue shaded regions do not describe our universe as the isocurvature density perturbations are too large.
know that ξ D is close to the virialization boundary, and from Figure 5 the form of the virialization boundary requires that p(f ) increases towards low f in the vicinity of the boundary. We assume, for simplicity, that the distribution can be approximated as a power law in the vicinity of the boundary, p(f ) = f n . We study which values of n lead to our proximity to the boundary and, in the Pre-Inflation case where θ is also scanning, we obtain the corresponding probability distributions for f that reproduce the observed DM abundance.
In order to obtain the posterior probability distribution for f we have to take into account measure and anthropic selection effects which modify the above simple power law behavior, and recall that our overall scheme includes the scanning of Λ. Following the discussion in section 3 that led to (16) and (17), after marginalizing over Λ the effective probability distribution for f and θ becomes dP ∝ θ(ξ D − ξ c ) where the virialization boundary is approximated by a θ function at ξ c = 0.5 ξ D0 . The integral (dθ) is present only for the Pre-Inflation case.
For the Post-Inflation cosmology, we parametrize where f 0 is the value for which the observed DM abundance is reproduced and has a theoretical uncertainty from the contribution of axion topological defects. Using the variablef the probability distribution reads wheref c 0.56 is the value off corresponding to ξ D = ξ c . In order to get a normalizable distribution for n > 1.2 we cut off the range off atf max = 2.15 × 10 3 , correspondent to ξ D = ξ max = 10 4 ξ D0 . This can be interpreted as an additional anthropic boundary at large ξ D related for instance to close stellar encounters [35,7]. The probability distributions for three different values of n are shown in Figure 7.
Turning to the Pre-Inflation scenario, the probability distribution is doubly differential, with also the initial misalignment angle scanning uniformly between −π and π. Hence dθ is included in (32) and the axion abundance in this case is given by (23). The 2-dimensional distribution in the (θ, f ) plane is not particularly illuminating. An important question to answer is about the support of the distribution in this case. We trade the variable f for ξ D , and as already discussed the DM abundance scans over the interval We require that we never scan over values of f greater than the Planck mass, which translates into the condition θ min = 10 −2 < θ < π , with θ min determined by ξ D (θ min , f = M Pl ) = ξ max . We use the probability distributions to compute the average value of ξ D as well as the 1σ confidence interval for both cosmologies, and we discuss the range of n that makes the observed abundance of DM in our universe typical. For Post-Inflation cosmology we change variable in the distribution (34) by using Eq. (33), and we derive a probability distribution for ξ D . For Pre-Inflation cosmology we start from the double differential distribution in Eq. (32), trade f for ξ D using Eq. (23) and then compute the average of ξ D in the domain described by Eqs. (35) and (36). The additional integration over θ in the Pre-Inflation case can be performed independently, thus in the two cases we have the same results, which are shown in Figure 8. We find that our universe is '1σ-typical' for −2.42 ≤ n ≤ 0.88. It is worth emphasizing that this range is not sensitive to the detailed choice of ξ max , and it stays unaffected for ξ max = (10 3 ÷ 10 6 ) ξ D0 .
In the Pre-Inflation case, for each value of n that successfully gives our proximity to the virialization boundary, we can obtain the probability distribution for measuring f in our universe given the observed value of ξ D0 . We thus proceed to integrate out θ from (32) imposing that the total axion density is the one observed today. The probability distribution for f thus reads Using (23) for ξ D and including anharmonic effects, numerical integration over θ yields the distributions of Figure 9, shown for three values of n. More generally, Figure 10 shows the average and 1σ ranges of f for all values of n that give our proximity to the virialization Only the smaller range of n above zero leads to the expectation of larger values of f . The case n = 0 is particularly interesting, occurring if the PQ scale is generated dynamically via a dimensional transmutation, and yields f ≈ (10 11 − 10 12 ) GeV. One interesting question to ask is how well an experiment like ADMX would perform in discovering the axion assuming it is the DM. The answer to this question depends on whether we live in the Pre-or Post-Inflation scenario and on the value of n. According to the experimental collaboration, the ADMX and ADMX-II experiments are going to be sensitive to the following ranges of f ADMX: 1.7 × 10 12 GeV f 3 × 10 12 GeV , ADMX-II: 3.4 × 10 11 GeV f 3 × 10 12 GeV , which are shown as shaded bands in Figures 9 and 10. Unfortunately neither of these two phases of the experiment is expected to cover a Post-Inflation axion, where 1.6 × 10 10 GeV f 1.6 × 10 11 GeV depending on whether the contri- bution from the decay of topological defects is taken into account or not. 8 The situation is different in the Pre-Inflation case, where a larger value of f can be accommodated by a small initial misalignment angle θ. However, the distribution for f is constrained by requiring our universe to be near the virialization boundary, and Figures 9 and 10 show that, for a wide range of allowed n, f is typically below the reach of these experiments. The probability for f to be in the (ADMX, ADMX-II) reach is shown in Figure 11(a,b). Thus for n not too far from zero there is hope that these experiments will make a positive discovery, but this should not detract from the intense need to design experiments to probe the lower f region.

The Thermal Axion Window
The observed DM abundance can be understood if the multiverse favors low values of f . If this is the case, we live close to the catastrophic boundary coming from requiring sufficient DM for density perturbations to go non-linear and halos to form by virialization. However, the argument presented in this Section is not quite complete because, for low enough f , sufficient axions are produced from thermal scattering for the axion DM density to rise above the virialization bound. This low f region is observationally excluded, for example from limits on axion emission from supernovae and from white dwarfs; we argue now that it is also anthropically disfavored.
The calculation of the density of thermal axions is performed in App. A. Skipping all the   Figure 11: Probability that f is in the region probed by ADMX (left) and ADMX-II (right) as a function of n. The gray region corresponds to the values of n for which our universe is not '1σ-typical'.
In order for this region not to be anthropically excluded from virialization, the axion density must above the catastrophic boundary ξ a > ξ c , so that f must be below a critical value The PQ breaking scale cannot be arbitrarily small, and the lower limit is set by a catastrophic boundary coming again from large scale structure. Small values of f make the axion strongly coupled to SM fields, and very short lived through the decay process a → γγ. If there is no DM around at recombination, baryons would not be able to fall into DM potential wells, and therefore all the perturbations below the Silk scale would be damped. This translates into the anthropic condition τ a→γγ ≥ t rec , with t rec the time when recombination happens. The axion life-time can be computed from the effective lagrangian in (21), and it results in The coupling c γ is model-dependent (see e.g. [44]), but typically of order one, so we set c γ = 1 in what follows. The recombination temperature is T rec 0.3 eV, for all purposes independent of the DM abundance. We use the Friedmann equation to obtain t rec ; for f ≤ f th c the axion is non-relativistic at recombination and the equivalence temperature satisfies T eq ≥ T rec , therefore we consider a matter dominated universe. We find that this imposes f > f th rec , where f th rec = 2.4 × 10 4 GeV .
There is indeed a thermal window where axion DM would not be excluded by virialization: f th rec ≤ f ≤ f th c . This low-f region is shown in Figure 12, where we plot both the axion relic density and the axion life time as a function of f (and the axion mass, upper axis). The shaded gray areas are excluded by virialization, but the region between them is in principle viable. Interestingly, the value of f such that the axion lifetime is of the order of the Hubble time falls into this range, f H 2.5 × 10 5 GeV.
There are anthropic reasons why this region is unlikely. Thermal axions produced for such low values of f are hot DM, and after they decouple they can free stream from overdense to underdense regions, damping primordial perturbations. The free streaming scale λ FS can be estimated by the expression (see e.g. [39]) λ FS 20 Mpc 10 eV m a .
This length has an associated mass, corresponding to the mass scale of the perturbation which are free streamed away, resulting in The mass density ρ FS has to account for the fact that axions can be unstable on a cosmological scale, thus we use the expression The final result is shown in Figure 13 and implies that, at the very least, large scale structure formation is very different from our universe. Throughout the thermal axion window perturbations on scales M < 10 9 M are destroyed by axion free streaming, so that population III stars fail to form. In about half of the window even galaxies of M = 10 12 M fail to form.
where we introduce the relativistic invariant phase space (49) depending on whether they are still in thermal equilibrium. The validity of this condition can be check in Figure 14a, where we plot the decoupling temperature as a function of f . So far we have not considered the possibility of an inflationary period wiping out the produced axions. However, even if this happens, axions can be generated again below the QCD phase transition from reactions with hadrons in the initial state. We consider the scattering process within the s-wave approximation, and upon using (48) we find the rate Γ a (ij ↔ ka) = n eq i n eq j n eq a |M ij→ka | 2 32πm i m j 1 − m k m i + m j 2 . (56) The dominant contributions to the total rate comes from the processes πη → πa , πK → Ka , πN → N a , ππ → πa .
The matrix elements M ij→ka for the first three processes can be computed by using the axion couplings to hadrons as derived in Refs. [42,43,44]. The process ππ → πa is not present if one considers the leading order Lagrangian, so we construct the next to leading order Lagrangian in the chiral expansion and compute the associated matrix element. Once we sum over all the contributions, we find that right below the QCD phase transition Γ a ≥ H for f ≤ 10 12 GeV, so that axions are kept in thermal equilibrium (or in the case of a low reheating temperature axions are regenerated). These reactions are effective until Γ a H, and we plot the decoupling temperature in Figure 14b. The decoupling temperature is not very sensitive to f , and has typical values of tens of MeV. We have included axion-hadron interactions that arise from the aGG coupling; the model-dependent couplings of axions to quarks could give order unity corrections to our results. To summarize, scattering processes above and below the QCD phase transition populate the universe with thermal axions. The reactions are effective until the decoupling temperatures shown in Figure 14b, and after that the axion number density just red shifts with the expansion. The number density today reads n a (T 0 ) = g * s (T 0 ) g * s (T dec ) The axion density is easily obtained as ρ a = m a n a , and we express it here as a fraction of the critical density Ω a h 2 = ρ a ρ cr h 2 0.12 3.4 × 10 5 GeV f = 0.12 m a 18 eV .