Superheavy Dark Matter from String Theory

Explicit string models which can realize inflation and low-energy supersymmetry are notoriously difficult to achieve. Given that sequestering requires very specific configurations, supersymmetric particles are in general expected to be very heavy implying that the neutralino dark matter should be overproduced in a standard thermal history. However, in this paper we point out that this is generically not the case since early matter domination driven by string moduli can dilute the dark matter abundance down to the observed value. We argue that generic features of string compactifications, namely a high supersymmetry breaking scale and late time epochs of modulus domination, might imply superheavy neutralino dark matter with mass around $10^{10}-10^{11}$ GeV. Interestingly, this is the right range to explain the recent detection of ultra-high-energy neutrinos by IceCube and ANITA via dark matter decay.


I. INTRODUCTION
While there are various lines of evidence for the existence of dark matter (DM) in the universe [1], the nature of DM remains a major problem at the interface of cosmology and particle physics. Weakly interacting massive particles (WIMPs) have long been a promising candidate and the focus of most direct, indirect and collider searches. In an attractive scenario, called the 'WIMP miracle', the DM relic abundance is obtained via thermal freeze-out in a radiation dominated (RD) universe for the nominal value of the DM annihilation rate σ ann v = 3 × 10 −26 cm 3 s −1 . This scenario, however, has been coming under increasing scrutiny by recent experiments, namely the Fermi-LAT results from observations of dwarf spheroidal galaxies [2] and newly discovered Milky Way satellites [3]. A recent analysis [4] has specifically ruled out thermal DM with a mass below 20 GeV in a model-independent way (unless there is P-wave annihilation or co-annihilation). Masses up to 100 GeV can be excluded if specific annihilation channels are considered.
The situation is different if the universe is not RD at the time of DM freeze-out [5]. This typically happens in non-standard thermal histories where the universe is not in a RD phase from inflationary reheating all the way to Big Bang nucleosynthesis (BBN) [6]. An important example is an epoch of early matter domination (EMD) driven by a component whose equation of state is the same as matter. This is a generic feature of early universe models arising from string theory constructions [7][8][9]. In this context, a string modulus is displaced from the minimum of its potential during inflation. Due to its long lifetime, the modulus dominates the energy density and gives rise to a period of EMD in the postinflationary history. The modulus eventually decays and a RD universe is established prior to BBN. Various production mechanisms during EMD can yield the correct DM abundance for both σ ann v < 3 × 10 −26 cm 3 s −1 and σ ann v > 3 × 10 −26 cm 3 s −1 [10].
Furthermore, the DM relic abundance can be completely decoupled from σ ann v if its main source is direct production from the decay of the component that drives an EMD phase [11]. In this scenario, the relic abundance depends on the branching fraction for decay to DM (hence the 'branching scenario' [12]) and the yield from the decay of the matter-like component. Non-thermal production of supersymmetric DM via the branching scenario has been studied in explicit string theory constructions where the volume modulus drives an epoch of EMD just before the onset of BBN [13]. A successful realization along these lines seems to be challenging for two reasons. First, the branching fraction of the volume modulus to DM is such that the correct abundance can be obtained for DM ∼ O(10) GeV. Second, the decay of the volume modulus typically produces dark radiation (DR) in addition to DM, and avoiding an excess of DR severely constrains the branching scenario [14].
However, in this work we shall show that the branching scenario could instead arise very generically in 4D string models with superheavy WIMPs. Several scenarios of supersymmetry breaking and inflation have already been realized in the context of string theory. Combining lowenergy supersymmetry with successful inflationary model building is notoriously hard to achieve [15]. The main reason is that the requirement of obtaining density perturbations of the correct size tends to fix the inflationary scale at relatively high energies. In turn, masses of the supersymmetric particles are also generically pushed to large values, typically at an intermediate scale around 10 10 − 10 11 GeV. A possible way to reconcile inflation with low-energy supersymmetry is to sequester the visible sector from the source of supersymmetry breaking in the bulk of the extra dimensions. Sequestered models, however, require a very specific brane configuration and Kähler metric for matter fields [16,17]. This solution therefore is not very generic. We note that statistical studies showed that high scale supersymmetry is a generic feature of the string landscape regardless of inflation [18,19].
Though the thermal DM scenario is known to overpro-duce superheavy WIMPs [20], the DM abundance may be diluted by epochs of EMD driven by string moduli. Hence two generic features of string compactifications, high-scale supersymmetry breaking and late time epochs of modulus domination, can successfully accommodate superheavy WIMPs with a mass around 10 10 −10 11 GeV.
Incidentally, if such a DM candidate is unstable and has the right coupling to neutrinos, its decay into very energetic neutrinos could provide a tantalizing explanation of the ultra-high-energy cosmic rays recently observed by IceCube and ANITA [21]. We will illustrate this general picture by presenting an explicit model that involves two periods of EMD. The first one is driven by inflaton oscillations at the end of which the inflaton mainly decays to DR in a hidden sector, and produces superheavy DM via its tiny coupling to the visible sector. A second stage of EMD is driven by the volume modulus, which is dominantly coupled to the visible sector and is lighter than the DM. As a result, this second EMD phase only dilutes the abundance of DM and DR produced in inflaton decay down to observationally acceptable values.
This paper is organized as follows. In Sec. II we briefly review the branching scenario for DM production. In Sec. III we discuss a successful framework for production of superheavy DM via the branching scenario. In Sec. IV we introduce an explicit string theory model for realizing this scenario. In Sec. V we identify the allowed parameter space of this model for a successful inflation and a correct DM abundance, and we present numerical results for the post-inflationary evolution for a benchmark point. We conclude in Sec. VI, and discuss generalized scenarios that involve more than one modulus in App. A.

II. BRANCHING SCENARIO: A BRIEF REVIEW
Let us consider a post-inflationary history that includes an EMD era driven by coherent oscillations or non-relativistic quanta of a long-lived scalar field ϕ with mass m ϕ and decay width Γ ϕ . The continuous decay of ϕ feeds radiation (assuming that decay products thermalize immediately) during the period that it dominates the energy density of the universe. The decay of ϕ completes when the Hubble expansion rate is H Γ ϕ , at which time the universe enters a RD phase. The resulting reheat temperature is T R = (90/π 2 g * ,R ) 1/4 (Γ ϕ M P ) 1/2 , where g * ,R is the number of relativistic degrees of freedom at T = T R .
The energy densities of ϕ and radiation, denoted by ρ ϕ and ρ R respectively, and the number density n χ of DM particles χ are found by solving the following system of Boltzmann equations: The first term on the right-hand side (RHS) of the last equation accounts for DM annihilation and inverse annihilation from the thermal bath ( σ ann v denotes the thermally-averaged annihilation/inverse annihilation rate). The second term accounts for direct production of DM from ϕ decay [22] (where Br χ is the number of DM particles produced per ϕ decay). Freeze-out/in of DM happens during the EMD epoch if T f > T R , where T f m χ /20 in the case of freeze-out and T f m χ /4 for freeze-in [23,24].
Assuming that freeze-out/in production is subdominant, the main contribution to the DM relic density comes from direct production at H Γ ϕ , and the number density of DM particles at this time is given by: The comoving number density of DM follows this expression, hence the name 'branching scenario' [11,12], provided that residual annihilation of DM particles to the thermal bath is inefficient. This will be the case if σ ann v n χ < Γ ϕ , where n χ is substituted from (2). Otherwise, partial annihilation will somewhat reduce the DM number density leading to the so-called 'annihilation scenario' of DM production [25,26]. The annihilation scenario can only be successful if σ ann v > 3 × 10 −26 cm 3 s −1 , which happens to be the case for weak-scale Wino and Higgsino DM. For small values of σ ann v , as in the case of Bino DM or for m χ 100 TeV, only the branching scenario can yield the correct DM abundance. After normalizing n χ in (2) by the entropy density s = 2π 2 g * ,R T 3 R /45 at T = T R , the DM relic abundance in the branching scenario is found to be: Here 3T R /4m ϕ is the yield factor that is related to dilution due to entropy released by ϕ decay. In order for the branching scenario to work, this must match the observed value: A natural question is if the branching scenario can be successfully realized in explicit particle physics models of the early universe. This issue has been discussed in the context of type IIB string compactifications where ϕ is the volume modulus [13]. In this case, we have T R /m ϕ (m ϕ /M P ) 1/2 . Also, for supersymmetric DM, three-body decays of ϕ result in a lower bound Br χ O(10 −3 ). Considering that T R 3 MeV (corresponding to m ϕ 50 TeV) is required for BBN, (3) and (4) imply that the correct DM abundance can be obtained for m χ O(10) GeV. Moreover, avoiding excessive production of DR, which typically accompanies DM production in string compactifications [27][28][29][30], seems to favor the annihilation scenario [14].

III. BRANCHING SCENARIO AND SUPERHEAVY DM
In this Section we lay down a framework for production of superheavy DM via the branching scenario. To overcome the challenges mentioned in Sec. II, we invoke two epochs of EMD driven by the inflaton and a modulus field respectively, as in generic string models. We also consider constraints from the cosmic microwave background (CMB) on such a scenario. In Sec. IV we shall present an explicit type IIB string model that successfully realizes this scenario (see also App. A for another explicit string model which realizes this scenario with an additional epoch of moduli domination).

A. Scenario with an epoch of modulus domination
The scenario we consider involves two periods of EMD driven by the inflaton σ and a modulus field φ in succession. Both of these fields behave as the field ϕ described in Sec. II. The inflaton σ is responsible for inflation at the end of which the Hubble expansion rate is H inf . The inflaton mass at the minimum of its potential is m σ and its couplings to the visible and hidden sectors are c vis /M P and c hid /M P respectively where c vis c hid . We will also assume that there is no stable non-relativistic particle in the hidden sector, so that the inflaton decay into hidden sector degrees of freedom just produces DR. Therefore, the inflaton decay rate into DR dominates over the one into visible sector particles since Γ σ→DR c 2 hid m 3 σ /M 2 P Γ σ→vis c 2 vis m 3 σ /M 2 P . We will also denote the total inflaton decay width as Γ σ = Γ σ→vis + Γ σ→DR .
The modulus φ has mass m φ < m σ . Its coupling to the visible sector is d vis /M P , while its coupling to the hidden sector is d hid /M P with d vis d hid . We will assume again that the modulus decay into the hidden sector produces just DR. This gives Γ φ→vis The total modulus decay width is instead Γ φ = Γ φ→vis +Γ φ→DR . We assume that m φ < m χ so that φ decay to DM is kinematically forbidden. The modulus acquires a displacement φ 0 from the minimum of its potential during inflation.
Below, we summarize the important stages of the postinflationary history in this scenario in chronological order: The universe is in an EMD phase driven by inflaton oscillations about the minimum of its potential. φ also starts oscillating at this stage and ρ φ = (φ 0 /M P ) 2 ρ σ . The inflaton decay completes at H Γ σ and mainly populates the hidden sector.
2-H D H < Γ σ : The universe is in a RD phase at this stage. The modulus oscillations behave like matter, and hence ρ φ is redshifted more slowly than ρ R . As a result, φ starts to dominate at H D (φ 0 /M P ) 4 Γ σ , which is the onset of a second phase of EMD.

3-Γ φ
H < H D : The universe is in a modulus-driven EMD epoch during this stage. The modulus decay completes when the Hubble expansion rate is H Γ φ and reheats the visible sector. This results in the formation of a RD universe prior to the onset of BBN.
The inflaton decay to the visible and hidden sectors produces DM and DR respectively. Given that m χ > m φ , the modulus decay dilutes both abundances and reproduces some amount of DR in the hidden sector. The number density of DM particles directly produced by the inflaton decay at H Γ φ is: where n σ = 3Γ 2 σ M 2 P /m σ is the inflaton number density at the end of stage 1, (a σ /a D ) 3 = (H D /Γ σ ) 3/2 is the number density redshift during stage 2, and (a D /a φ ) 3 = (Γ φ /H D ) 2 is the number density redshift during stage 3.
If DM is the lightest R-parity odd particle in the visible sector, we have: The first factor on the RHS of this expression is the fraction of σ quanta that decay to the visible sector. The second factor is the ratio of the number of R-parity odd particles (which subsequently decay to DM) to the total number of particles in the visible sector produced per σ decay. In the explicit example that we discuss later, the σ decay into the visible sector mainly occurs through two-body decays to gauge fields. Two-body decays to Rparity odd particles are highly suppressed, but they are produced via three-body decays including one gauge field and two gauginos resulting in Br vis,odd 10 −3 (which is essentially a phase space factor) [12]. Therefore, after normalizing n χ by the entropy density s, we find: where: with g * ,R denoting the number of relativistic degrees of freedom in the visible sector at T = T R , and Y φ ≡ φ 0 /M P . Regarding DR, its energy density at H Γ φ is: where ρ σ 3Γ 2 σ M 2 P , (a σ /a D ) 4 = (H D /Γ σ ) 2 is the energy density redshift during stage 2, (a D /a φ ) 4 = (Γ φ /H D ) 8/3 is the energy density redshift during stage 3, and ρ φ 3Γ 2 φ M 2 P . Hence, the final fractional energy density of DR is given by: This ratio must be small enough to satisfy the observational constraints on the DR abundance.

B. Constraints from CMB
Inflation is the dominant paradigm for generating the almost scale-invariant perturbations. The number of efoldings between the time when perturbations of a given wavelength exit the horizon and the end of inflation depends on the scale of inflation as well as the postinflationary thermal history. One or more periods of EMD change the number of e-foldings from that in a standard thermal history.
In the scenario discussed in Sec. III A, the number of e-foldings of inflation between the time when the pivot scale k * = 0.05 Mpc −1 left the horizon and the end of inflation can be written as [31,32]: where r is the tensor-to-scalar ratio and: Here N reh and N φ denote the duration of EMD phases from inflationary reheating and modulus domination (stages 1 and 3 above) respectively. This results in: In important universality classes of inflation, the scalar spectral index n s is related to N e through a simple relation [33]: For example, in the Starobinsky model and Higgs inflation, as well as the specific model of string inflation that we will discuss later, a = 2. This then leads to: This implies that: where n s,min is the minimum value in the 2σ region allowed by Planck data [34]. For a given model of inflation where H inf is known, this in turn sets an upper bound on Y 4 φ Γ −1 φ through (13). On the other hand, for known inflaton parameters m σ and Γ σ , (7) and (10) result in a lower bound on Y 4 φ Γ −1 φ in order not to overproduce DM and DR in our scenario. Therefore, obtaining the correct abundance of DM (while avoiding an excessive production of DR) and getting an acceptable value of n s constrain the epoch of modulus domination in opposite ways. 1 This can be understood intuitively as follows. While diluting the abundance of DM and DR produced from inflaton decay to acceptable levels requires a long enough bout of modulus domination, satisfying the lower bound on n s limits the duration of that period from above.

IV. A STRING MODEL WITH AN EPOCH OF MODULUS DOMINATION
In this Section we shall present an explicit string model which successfully realizes inflation and superheavy DM via the branching scenario with an epoch of modulus domination.

A. The setup
We consider a type IIB model with 3 Kähler moduli T i = τ i + ic i , i = 1, ..., 3 and a Calabi-Yau volume of the form: The visible sector is realized via a stack of D7-branes wrapped around the 4-cycle whose volume is controlled by τ vis , while inflation is driven by the modulus τ inf as in Kähler moduli inflation [36]. A hidden sector lives instead on a stack of D7-branes wrapped around the 4cycle whose volume is given by τ inf . The structure of the effective supergravity theory is determined by the Kähler potential K and the superpotential W . K is given by: where g s is the string coupling and ξ is an O(1) coefficient which controls α corrections [37] beyond the tree-level expression. W instead reads: where W 0 ∼ O(10 − 100) is the tree-level contribution, while the terms proportional to A vis and A inf are nonperturbative effects [38] (all A's and a's are expected to be O(1) constants). Moduli stabilization produces a typical LVS minimum [39] at exponentially large volume in string units, V τ 3/2 big ∼ e 1/gs , while the two blow-up modes are fixed at smaller values τ vis ∼ τ inf ∼ 1/g s ∼ O(10), where we take the string coupling in the perturbative regime g s 0.1. Notice that τ vis sets the value of the visible sector gauge coupling α −1 vis = 4πg −2 vis = τ vis ∼ O(10) which turns out to be in the appropriate phenomenological regime.
Moduli stabilization proceeds as follows: at leading order in a 1/V expansion, non-perturbative corrections to W combined with α corrections to K stabilize V, τ vis , c vis , τ inf and c inf , leaving 1 flat direction parameterized by the axion c big . 2 This axion turns out to be ultra-light since it receives a tiny mass due to additional T big -dependent non-perturbative corrections to W . Thus c big plays the role of hidden sector dark radiation. This system admits a non-supersymmetric AdS minimum which can however be uplifted to dS via several possible mechanisms (anti D3-branes [42], T-branes [43], non-perturbative effects at singularities [44], non-zero Fterms of the complex structure moduli [45]).

B. Moduli mass spectrum
The determination of the moduli mass spectrum and couplings to both visible and hidden sector fields requires first to go to canonically normalized fields. Following the notation of Sec. III, we will denote them as: (i) σ for τ inf since this modulus plays the role of the inflaton; (ii) φ for τ big since this modulus will give rise to an EMD epoch after the end of inflation; and (iii) a DR for the closed string axion c big which behaves as dark radiation. Defining: the mass spectrum of the relevant moduli around the minimum becomes [46,47] (see [48] for the correct normalization factor κ): This setup allows to realize Kähler moduli inflation [36] where the inflaton is σ since this modulus becomes much lighter than H m φ as soon as it is displaced from its minimum. τ vis , c vis , and c inf are heavy spectator fields which do not get displaced during inflation since their mass is of the same order of the mass of σ around the minimum, and so it is much larger than H. On the other hand, all the other moduli get displaced from their minimum during inflation. We shall focus just on φ since the axion a DR remains almost massless and behaves as a source of extra dark radiation. We shall also denote the displacement of the canonically normalized light Kähler [9]. Due to this displacement during inflation, φ gives rise to a period of modulus domination. Moreover supersymmetry is broken due to non-zero F-terms of the Kähler moduli which generate a gravitino mass m 3/2 together with gaugino and scalar masses of order [49]: Taking the DM mass of the same order as the soft terms, m χ m 0 M 1/2 , we realize that: which ensures that DM cannot be reproduced from the decay of the light modulus φ.
Notice that in order to avoid any cosmological moduli problem, the mass of φ has to be m φ O(50) TeV. Using (21) and setting g s 0.1 and 1 W 0 100, this gives the bound 5 × 10 −9 − 10 −8 1 which, when translated in terms of the overall volume, becomes 1 V 10 8 −5×10 9 . This, in turn, produces a scenario of superheavy DM since it sets a lower bound on the DM mass of order m χ 10 10 − 10 11 GeV. As we shall see in the Sec. V, values of V below 10 8 − 5 × 10 9 are also required to generate, during inflation, the observed value of the amplitude of the density perturbations.

C. Hidden sector configuration
Let us comment a bit more on the configuration of the hidden sector D7-stack wrapping τ inf . This has to provide a non-perturbative contribution to the superpotential which generates the inflationary potential, and be such that the inflaton decay into the hidden sector produces just relativistic degrees of freedom without additional contributions to the DM abundance. This dark radiation component is subsequently diluted by the decay of the lightest modulus. If the hidden sector is a supersymmetric SU (N c ) theory with N f flavors, it would confine if N f < N c . The corresponding scale of strong dynamics Λ can be shown to be above the inflaton mass, m σ < Λ [46], and so σ cannot decay into glueballs (gg), 'gluinoballs' (gg), and 'glueballinos' (gg) since they all develop a mass of order Λ. Hence we need to discard the pure SYM case. For N f > 0 with soft supersymmetry breaking terms, squarks and quarks form scalar and fermionic condensates which all develop a mass of order m 0 M 1/2 m σ [50], except for pion-like mesons which are exactly massless in the absence of a supersymmetric quark mass term in W . Therefore σ could decay into these heavy condensates but some of them would be stable in the absence of EW interactions. We conclude that the hidden sector cannot be a simple SU (N c ) theory with N f flavors. The best configuration for the hidden sector is instead a copy of the visible sector, i.e. an MSSM-like hidden sector, with however 3 differences with respect to the visible sector: (i) the scale of strong dynamics Λ is much higher than in ordinary QCD; (ii) R-parity is not conserved so that hidden protons are unstable; (iii) the mass of the hidden electrons is very small so that they are still relativistic, like neutrinos. In this scenario, all hidden degrees of freedom produced from the inflaton decay eventually decay into hidden massless gauge bosons or hidden relativistic matter fermions.
A final requirement is the absence of any leakage of energy between hidden and visible sector degrees of freedom due to kinetic mixing between U (1)s or a possible moduli portal. The first option can be avoided by construction if the hidden gauge group does not contain any Abelian U (1) factor. 3 On the other hand, a moduli portal between the two sectors could be created by the volume modulus φ. However, we expect this effect to be negligible since, as we shall see in Sec. IV D, this field couples only with Planckian strength to both sectors, and so any leakage would be proportional to (1/M P ) 4 .

D. Moduli couplings and decay rates
Due to the geometric separation in the extra dimensions between τ vis (which supports the visible sector D7 stack) and τ inf (which supports a hidden sector D7 stack), the coupling of the canonically normalized inflaton σ to hidden sector gauge bosons is much stronger than the one to visible sector gauge fields [46]: with: Notice that the interactions in (24) provide the main contributions to the inflaton decay rate to both visible and hidden degrees of freedom. In fact, since m 0 M 1/2 m σ , the inflaton decay into supersymmetric partners is mass suppressed. The same consideration applies to the inflaton decay into both visible and hidden sector matter fermions. The decay rate into Higgses is also mass suppressed except for the case of a Giudice-Masiero interaction in K which we assume to be absent. 4 This implies the following important relation for the determination of the DM abundance using (7): where we included also the number of visible and hidden sector gauge bosons denoted respectively as N g and N hid g . For an MSSM-like visible sector we have N g = 12 while N hid g is a model-dependent parameter which can also be larger than N g .
On the other hand the light modulus φ can decay to: • Hidden sector gauge bosons: • Visible sector gauge bosons: • Visible sector Higgs h 0 and would-be Goldstone bosons G 0 and G ± [29]: For 1, the decay rate of φ into both hidden and visible gauge bosons is suppressed. On the other hand, the decay of φ into ultra-light bulk axions could give rise to extra dark radiation which needs to be in agreement with present observational bounds [52]. This sets a lower bound on Z of order (neglecting DR from inflaton decay since this is diluted by the decay of φ) [29]: Therefore the ratio Γ φ /Γ φ→vis which appears in (7) for the final DM abundance looks like: and the reheating temperature T R in (8) can be derived from the following decay width: Finally, the remaining quantities which are crucial to derive the fractional energy density of DR using (10) are: and: E. Consistency of the branching scenario In this paper we are considering a branching scenario for DM production from inflaton decay. This is generically the case for superheavy WIMP DM since the corresponding annihilation rate would be too small to realize the so-called non-thermal annihilation scenario. However the computation of the DM relic density relies on the assumption that the standard thermal freeze-out mechanism cannot occur. This is true if the visible sector reheating temperature after the inflaton decay T vis R,inf is below T f , where T f m χ /20 in the case of freeze-out (and T f m χ /4 for freeze-in). We shall now show that this is indeed the case in our model.
In standard supersymmetric scenarios, the LSP mass is expected to be of order the soft mass. As we have already seen, the DM mass is therefore slightly below the inflaton mass, m χ m σ /(ln ) 2 < m σ . Notice that this feature is not a peculiarity of our model but it is a generic characteristic of string compactifications since, whenever the 4-cycle supporting the visible sector is stabilized in the geometric regime, the visible sector is always not sequestered from the source of supersymmetry breaking in the bulk. Thus the soft terms, and the DM mass, turn out to be of the same order as the gravitino mass, which sets also the order of magnitude of the mass of generic moduli (up to possible | ln | suppression factors).
Therefore we shall consider T f m σ /[20 (ln ) 2 ]. On the other hand, the visible sector reheating temperature reads: where we used c hid c −1 vis , and N hid g N g = 12 .
Hence we obtain T vis R,inf < T f provided that: This is indeed the case for 1 and κ 1, which guarantees the consistency of the branching scenario. This will be confirmed in Sec. V B which presents a numerical analysis of the cosmological evolution of our model.

V. COSMOLOGY OF THE STRING MODEL
In this Section we shall first determine the values of the microscopic parameters which give the right amplitude of the density perturbations and the correct DM abundance, finding a DM mass around 10 10 -10 11 GeV. We shall then perform a numerical analysis of the cosmological evolution of our string model with an epoch of modulus domination.

A. Inflationary observables and DM abundance
Let us derive the allowed DM mass range in a single modulus cosmology. We achieve a rather precise prediction by imposing a combination of observational and geometrical constraints. We start with the expression for the number of e-foldings between horizon exit and the end of inflation [9]: Here r is the tensor-to-scalar ratio, N reh is the duration of the reheating period due to the inflaton σ, and N φ is the duration of the EMD epoch due to the modulus φ. Note that we have set the equation-of-state parameter w equal to zero during inflationary reheating. Also ρ σ,start is the energy density at horizon exit, while ρ σ,end is the energy density at the end of inflation. Let us rewrite (34) in terms of fundamental parameters. The duration of the reheating period is: where H σ,end is the Hubble rate at the end of inflation, which is given by [9]: Combining (36) with the inflaton decay rate (30) gives: The duration of modulus domination is given by: Following again [9], the tensor-to-scalar ratio can be expressed as: Noting that the amplitude of the density perturbations can be written as A s = 2 3π 2 r ρσ,start M 4 P , we get: where we have set Z = 2 and used ln 10 10 A s = 3.044 [34]. To proceed further, we need the relation between the inflaton τ inf , the volume V, and the number of efoldings N e that matches the observed value of A s . This reads [9]: Given that τ inf describes the volume of a local 4-cycle, we have to impose the geometrical constraint V 2/3 τ b τ inf which guarantees that the effective field theory is under control. We can implement this constraint as V 2/3 λ τ inf , where λ 1 is a tunable parameter that determines the hierarchy between the overall volume V and the volume of the blow-up mode τ inf . This gives us the final expression: with α = 1.15 × 10 −11 1 2π and N e given in (40). Let us briefly summarize the procedure that we shall follow to derive the DM mass corresponding to the observed DM abundance: • We extract from (42) W 0 as a function of V. This step encodes in W 0 (V) the information of the amplitude of the density perturbations and the geometrical relation between V 2/3 and τ inf . We also set a natural bound on W 0 by constraining it to be in the range O(1 − 10 3 ).
• We extract the value of V by matching the expressions for the observed and predicted DM abundances. We perform this step for each of the 1444 initial parameter combinations.
• We compute the DM mass for those parameter combinations that allow for the correct DM abundance. In Fig. 1 we present all data points in the (W 0 , V) plane which reproduce the observed amplitude of the density perturbations, respect our geometrical constraints, and yield the correct DM abundance. Approximately 72% of our initial parameter space leads to a consistent solution. Notice that, although each point corresponds to different values of g s , Y φ , and λ, the resulting DM mass is always in the range 10 10 -10 11 GeV, giving a robust prediction which is almost independent of the variation of the underlying parameters.
The accumulation around the origin and the jet-like structures in the distribution of the data points can be understood from Fig. 2 where we split the points shown in Fig. 1 into two sets with, respectively, g s = 0.001 − 0.009 and g s = 0.01 − 0.1. Moreover black dots correspond to λ = 10 4 , red to λ = 10 3 , blue to λ = 10 2 , and green to λ = 10. The plot for smaller values of g s shows clearly that the four jet structures correspond to different values of λ. This behavior is a direct consequence of (42) which implies that λ determines the slope of the function W 0 (V). On the other hand, the plot for larger values of g s features a larger density at smaller values of W 0 and V. This behavior is a consequence of the consistency of the branching scenario. In fact, in order for (26) to hold, smaller values of the volume must lead to an increase in g s . It is worth mentioning also that around 71% of the data points correspond to g s ∈ [0.01, 0.1], whereas 29% of the acceptable parameter space correspond to g s ∈ [10 −3 , 0.01]. In Fig. 3 we present a similar analysis, this time splitting all data points from Fig. 1 into two sets, depending on the value of the misalignment Y φ . We observe again the same jet structure depending on the value of the parameter λ. An important observation here is the slight rotation of the data-point cone towards the W 0 axis if we increase Y φ . This behavior is mainly driven by the DM abundance constraint formulated in (7). The abundance scales like ∼ Y φ −2 V −13/4 . Hence, in order to match the right abundance, a smaller volume must be compensated by a larger misalignment Y φ .
Understanding the behavior of our data set as a function of the underlying parameters is important in order to understand the distribution of the scalar spectral index n s . For each point in Fig. 1, we calculated the resulting value of n s . All obtained values are within the 2-and 3σ range [34], as can be seen from Fig. 4.
In Fig. 4 each black dot corresponds to a scalar spectral index within the 2σ range, i.e. 0.9565 < n s < 0.9733, Around 42% of the data points are in the lower Y φ regime, while around 58% are in the upper regime. Black points correspond to λ = 10 4 , red to λ = 10 3 , blue to λ = 10 2 , and green to λ = 10. while each red dot has n s in the 3σ range, i.e. 0.9523 < n s < 0.9775, but outside 2σ. By comparing the two distributions, we see that the data points corresponding to 2σ seem to be accumulating at the origin and their cone is slightly rotated towards the volume axis. From this observation we conclude that phenomenologically more acceptable values of n s drive the string coupling g s to larger values and the misalignment Y φ to smaller ones. A crucial observation is that the parameter values in our data set naturally respect the relation between the volume and the underlying model parameters at the minimum of the scalar potential [39]: , for natural O(1) values of the microscopic parameters a inf , A inf and ξ.
In Sec. V B we shall perform a more in-depth numerical analysis of the cosmological evolution, using the benchmark parameters listed in Tab. I.

B. Numerical analysis of cosmological evolution
We perform a numerical analysis of the cosmological evolution of our scenario by solving the coupled set of Boltzmann equations for the various cosmological components (see App. A for a scenario with two moduli). We begin the numerical evolution at H H inf , with both σ and φ oscillating, and other components highly subdominant. The Boltzmann equations for our single-modulus scenario are as follows: where the Hubble rate H is given by the sum of all energy density components, and the various decay rates are given in (26), (29) and (30) using the benchmark values of Tab. I. σ ann v denotes the thermally averaged rate for χ production from/annihilation to the thermal bath with the average energy per χ particle approximated as E χ ≈ m 2 χ + 9T 2 vis [23]. Here, we take σ ann v ≈ α 2 χ /m 2 χ with α χ ∼ 0.1. This happens to be the case, for example, for Higgsino and Wino DM [53]. However, because thermal production is subdominant in our scenarios, the exact form of σ ann v , including possible temperature dependence, is not really important. For typical DM masses in our scenarios, m χ ∼ 10 10 -10 11 GeV, we obtain values of σ ann v in the freeze-in regime. Finally, the DM equilibrium number density, relevant for thermal production, is given by: A sample numerical solution of (43)- (47) is shown in Fig. 5. As the evolution proceeds, DM, dark radiation, and ordinary radiation are continually produced by inflaton decay until H Γ σ , at which point inflaton decay completes. This begins an era of hidden-radiation domination which lasts until the light modulus φ overcomes the energy density of hidden radiation. From here until the time when H Γ φ , we have a period of EMD driven by the modulus, which is then followed by the standard period of radiation domination once the modulus decay completes. The DM abundance is set by the inflaton decay at H Γ σ and simply redshifts through the remaining cosmological history. For typical values of the parameters in our scenario, the maximum visible sector temperature established during inflationary reheating is smaller than the DM mass, such that thermal production of DM occurs on the Boltzmann tail of the equilibrium distribution, rendering the thermal contribution irrelevant. 5 Fig.  6 shows the visible sector temperature as a function of the scale factor for the cosmological history shown in Fig.  5, where we have assumed a smooth function for the temperature dependence of the relativistic degrees of freedom in the visible sector.
One comment is in order at this point. Our calculation of freeze-in production of DM in (47) assumes instantaneous thermalization of inflaton decay products in the visible sector. In fact, it holds as long as the visible sector reaches thermal equilibrium at a temperature T > T f . However, due to the small number density of inflaton decay products in the visible sector, thermalization may be significantly delayed (for example, see [54,55]).
If T < T f at the time of thermalization, then thermal production of DM will be completely negligible. Before thermal equilibrium is established, DM production from inflaton decay products is kinematically possible due to their typical mass hierarchy m σ m χ [56][57][58][59]. However, by conservation of energy, the number density of these particles is much smaller than what it would be in thermal equilibrium. We have checked that DM production during thermalization is a few orders of magnitude smaller than that from direct inflaton decay, for the parameters shown in Tab. I, even if the visible sector is not thermalized until H Γ σ .

VI. CONCLUSIONS
In this paper we have argued that two generic features of string compactifications, a high supersymmetry breaking scale (which is favored by both statistical arguments [18,19] and by the requirement of a viable inflationary model building [15]) and the presence of light moduli which drive epochs of EMD [7][8][9], lead typically to superheavy WIMP DM with mass around the intermediate scale. This scenario has not received significant attention so far because DM with a mass in the 10 10 -10 11 GeV range is inevitably overproduced in a standard thermal history. However, if DM is produced non-thermally from the decay of the inflaton and it is subsequently diluted by the decay of long-lived string moduli which are so light that their decay does not reproduce DM, one can obtain the observed abundance in this so-called branching scenario even for very large DM masses. This does not just account for the non-observation of supersymmetry and WIMPs at colliders, but it may also provide a natural explanation of the origin of ultra-high-energy cosmic rays recently observed by IceCube and ANITA, if DM is unstable and has the right coupling to neutrinos [21].
We illustrated this general picture by presenting two explicit 4D string models which lead to superheavy WIMP DM. The first model is described in Sec. IV and features a single epoch of modulus domination, while App. A gives all the details of a different model with two epochs of EMD driven by two different light moduli. It turns out that in both cases the observed DM abundance can be obtained for a mass around 10 10 -10 11 GeV. The main virtue of both models is the possibility to follow their entire cosmological evolution from inflation to the final reheating (due to the decay of the lightest modulus) that establishes a RD universe before the onset of BBN. This can be achieved by focusing on type IIB LVS string models where the exponentially large volume of the extra dimensions allows to keep control over the 4D low-energy effective field theory. Hence all moduli masses and couplings to both visible and hidden sector degrees of freedom can be computed in detail. Moreover one can build 4D models which can realize inflation, supersymmetry breaking and a chiral MSSM-like visible sector on D-branes (see [61][62][63][64][65][66] for explicit Calabi-Yau models with all these features).
We followed the entire cosmological evolution in both models using analytical and numerical tools. This allowed us to combine various constraints coming from both theoretical and phenomenological considerations. Interestingly, we derived the ranges of the microscopic parameters in a regime where geometrical constraints on the underlying extra-dimensional construction are respected, which yield the observed DM abundance as well as the correct value of inflationary observables, namely the amplitude of the density perturbations and the scalar spectral index.
Future investigations could include more formal aspects as well as more phenomenological implications of our findings. From the formal point of view, it would be very interesting to investigate how generic superheavy WIMP DM is from the string landscape point of view, for example comparing this scenario to the case of fuzzy DM [67] which has also been claimed to be a natural outcome of string models due to the ubiquitous presence of ultra-light axions [68][69][70][71]. On the other hand, from the phenomenological side, it is crucial to understand how a superheavy DM could be detected in actual observations, for example establishing in more detail the possible connection of our results with the production of very energetic cosmic rays from DM decay. We leave all these intriguing possibilities for future work. String compactifications are in general characterized by a large number of moduli and by a leading order noscale structure that makes some of them (in general more than one) lighter than expected, i.e. m φ m 3/2 [72]. Hence one might expect to have several epochs of modulus domination in the post-inflationary history. One could have either multiple eras of EMD separated by intermediate phases of RD (if there is a hierarchy among the initial displacements of the moduli), or a single extended EMD epoch (if the initial displacements of all moduli are of the same order). In what follows we shall investigate the features of this general scenario by focusing on the illustrative example with two light moduli.

Branching scenario with two epochs of EMD
Let us present a branching scenario that involves two moduli φ 1 and φ 2 . The two moduli have masses m φ1 and m φ2 with m φ2 m φ1 . They mainly decay to the visible sector with respective couplings c 1 /M P and c 2 /M P , leading to decay widths Γ φ1→vis c 2 1 m 3 φ1 /M 2 P and Γ φ2→vis c 2 2 m 3 φ2 /M 2 P (we assume that their decays to hidden sector particles are suppressed and produce just a small amount of DR). Assuming that m φ1 < m χ , the decay of φ 1 and φ 2 to DM will be kinematically forbidden. Therefore their decay only dilutes the abundance of DM and DR produced from the inflaton decay. The initial displacements of φ 1 and φ 2 from the minimum of their potential are φ 1,0 and φ 2,0 respectively. We assume φ 1,0 φ 2,0 so that each modulus can dominate the energy density of the universe for a period of time.
The important stages of the post-inflationary history in this scenario (in chronological order) are as follows: The universe is in an EMD phase driven by inflaton oscillations about the minimum of its potential. Both moduli also start oscillating at this stage with respective energy densities ρ φ1 = (φ 1,0 /M P ) 2 ρ σ and ρ φ2 = (φ 2,0 /M P ) 2 ρ σ . The inflaton decay completes at H Γ σ mainly populating the hidden sector.
2-H D,1 H < Γ σ : The universe is in a RD phase at this stage. Moduli oscillations behave like matter, and hence ρ φ1,2 redshifted more slowly than ρ R . Since ρ φ1 > ρ φ2 , φ 1 starts to dominate at H D,1 (φ 1,0 /M P ) 4 Γ σ , which is the onset of a second phase of EMD.

3-Γ φ1
H < H D,1 : The universe is in an EMD epoch that is driven by φ 1 during this stage. Decay of φ 1 completes when the Hubble expansion rate is H Γ φ1 (Γ φ1 denotes the total decay rate of φ 1 ) and reheats the visible sector. This leads to the formation of a RD universe.
4-H D,2 H < Γ φ1 : The universe is in an intermediate RD phase during this stage. Since ρ φ2 is redshifted more slowly than ρ R , φ 2 starts to dominate when the Hubble expansion rate is H D,2 (φ 2,0 /φ 1,0 ) 4 Γ φ1 , which is the onset of another epoch of EMD. 5-Γ φ2 H < H D,2 : The universe is in a third phase of EMD that is driven by φ 2 . The decay of φ 2 completes when the Hubble expansion rate is H Γ φ2 (Γ φ2 is the total decay rate of φ 2 ), at which time the universe enters the final RD phase prior to the onset of BBN.
One point to note is that H D,2 Γ φ1 if φ 2,0 φ 1,0 . In this case, φ 2 dominates the energy density of the universe as soon as the decay of φ 1 completes. Stage 4 above thus effectively disappears and there is a direct transition from the first EMD era (stage 3) to the second one (stage 5), implying a single extended epoch of EMD driven by the two moduli φ 1 and φ 2 . 6 Let us now estimate the relic abundance of DM in this scenario. The number density of DM particles at H Γ φ2 is given by: (A1) This is similar to (5) for the case with a single epoch of modulus domination. The last two terms on the RHS, which are new, account for the dilution of the number density in stages 4 and 5 above respectively. After using the scaling of a with time in stages 4 and 5 above, and normalizing n χ by the entropy density s at H Γ φ2 , we find: where: T R,2 = 90 π 2 g * ,R,2 with g * ,R,2 denoting the number of relativistic degrees of freedom in the visible sector at T = T R,2 , and Y φ2 ≡ φ 2,0 /M P . The energy density of DR at H Γ φ2 is given by (assuming that no DR is produced from φ 1 decay): This is similar to (9) with two additional terms on the RHS that account for the energy density redshift in stages 4 and 5 respectively. Thus the fractional energy density of DR is given by: (A5) Some comments are in order at this point. It is seen from (A2) and (A5) that final abundance of DM and DR in the two modulus scenario depends only on the initial amplitude and decay width of the second modulus φ 2 . This can be understood as follows. For fixed Y φ2 , varying Y φ1 ≡ φ 1,0 /M P (as long as Y φ1 Y φ2 ) affects the two epochs of modulus domination in opposite ways. Increasing (decreasing) Y φ1 makes stage 4 longer (shorter) and stage 3 shorter (longer) by the same factor. A similar thing happens by decreasing (increasing) Γ φ1 with Γ φ2 fixed (as long as Γ φ1 Γ φ2 ). As a result, the combined dilution factor from two epochs of modulus domination does not depend on the parameters of φ 1 .
That said, it is helpful to compare the DM and DR abundance with the previous scenario when there is one epoch of modulus domination. We can rewrite (A2) in terms of (7) as follows (in the limit where the production of DR from the decay of the lightest modulus is completely negligible): Similarly, (A5) can be written in terms of (10): It is seen from (A6) and (A7) that the maximum dilution in the scenario with two moduli is achieved for Y φ2 Y φ1 and m φ2 m φ1 . 7 As pointed out earlier, in this case there is a single extended epoch of EMD consisting of stages 3 and 5 that are not separated by an intermediate RD phase.

A string model with two epochs of modulus domination a. The setup
We now focus on a type IIB model which can allow for two epochs of modulus domination. This model shares the same features with the model discussed in Sec. IV but it also gives rise to novel phenomenological properties. The Calabi-Yau volume now takes the form: where again τ inf drives inflation and it is wrapped by a hidden sector D7 stack as in the model presented in Sec. IV. However now the second blow-up mode, here denoted as τ np , is just responsible for generating non-perturbative effects needed for moduli stabilization but it does not support the visible sector stack of D7 branes. 8 In fact, in this model the requirement to avoid dark radiation overproduction from the decay of the lightest modulus forces the visible D7 stack to be wrapped around the 4-cycle whose volume is controlled by τ vis [30].
In this case the 4D low-energy supergravity theory is determined by the following Kähler potential: where K gs denotes string loop corrections [73][74][75] which have been shown to be V-suppressed with respect to the leading α correction proportional to ξ [37]. The superpotential instead looks like: As in Sec. IV, at leading order in 1/V 1, nonperturbative corrections to W combined with α corrections to K produce an LVS minimum with 5 stabilized moduli: V √ τ vis τ big ∼ e 1/gs , τ np ∼ τ inf ∼ 1/g s ∼ O(10) together with the 2 corresponding axions c np and c inf . However at this level of approximation there are still 3 flat directions which can be parameterized by τ vis , c vis , and c big . The visible sector modulus τ vis is fixed at subleading order by the string loop contribution to the Kähler potential K gs at [76]: where λ loop is a tunable combination of the coefficients of g s corrections to K. Notice that the requirement of reproducing the observed value of the visible sector gauge coupling, α −1 vis = 4πg −2 vis = τ vis ∼ O(10 − 100), leads necessarily to an anisotropic shape of the extra dimensions since the exponentially large Calabi-Yau volume V √ τ vis τ big is now controlled by 2 cycles but with τ big ∼ e 1/gs τ vis ∼ 1/g s . Finally, the two axions c vis and c big receive tiny masses due to additional T vis -and T big -dependent non-perturbative corrections to W . Thus both c vis and c big are in general ultra-light and play the role of hidden sector dark radiation.

b. Moduli mass spectrum
The mass spectrum of the relevant moduli around the minimum becomes: for , g s 1 where σ, φ 1 , φ 2 , a DR1 , and a DR1 are the canonically normalized fields corresponding respectively to τ inf , τ big , τ vis , c big , and c vis . As in Sec. IV, σ plays the role of the inflaton. This field, when displaced from its minimum, becomes exponentially lighter than the Hubble constant during inflation which is set by the mass of φ 1 : H m φ1 . The 3 fields τ np , c inf , and c np are instead heavy spectator fields, while φ 1 and φ 2 get displaced from their minimum during inflation, and so give rise to 2 epochs of EMD. On the other hand, the 2 ultra-light axions a DR1 and a DR2 yield extra contributions to N eff . Moreover the gravitino mass and the soft terms are still given by (22). Hence for m χ m 0 M 1/2 , we conclude that DM cannot be reproduced from the decay of any of the 2 light moduli since: Requiring m φ2 O(50) TeV in order to avoid any cosmological moduli problem together with τ vis ∼ O(100) in order to reproduce the observed value of the visible sector gauge coupling, corresponds to 1 V 5 × 10 7 − 10 9 for g s 0.1 and 1 W 0 100. Therefore DM is necessarily superheavy since m χ 10 11 GeV. Notice that values of the overall volume below 5 × 10 7 − 10 9 are also required to match inflationary observables like the amplitude of primordial fluctuations [36].

c. Moduli couplings and decay rates
The configuration of the hidden sector D7-stack wrapped around τ inf is the same as the one described in Sec. IV C. Moreover, also the couplings of the inflaton σ to hidden and visible gauge bosons are still given by (25). Hence the ratio Γ σ→vis /Γ σ is also still given by (26) and the inflaton decay width Γ σ takes the same form as (30).
On the other hand the light modulus φ 1 decays mainly into visible sector gauge bosons with decay rate [30,46]: (A14) where γ ≥ 1 is a microscopic parameter which depends on the gauge flux on the visible sector D7-stack (in particular γ = 1 for a fluxless D7-stack while γ > 1 for non-zero gauge fluxes) [30]. The decay of φ 1 produces also axionic dark radiation which is however suppressed for γ > 1 and gets diluted by the decay of φ 2 . In what follows we shall therefore neglect Γ φ1→DR .
The final modulus to decay is φ 2 whose main decay channels are [30]: • Dark radiation bulk axions: • Visible sector gauge bosons: where we have set again N g = 12.
The amount of axionic dark radiation produced from φ 2 decay is controlled by the underlying parameter γ: showing how for γ ≥ 1 this model naturally satisfies present observational bounds since it yields ∆N eff 0.6. The relevant quantities to compute the final DM abundance using (A2) and (A3) can be derived from the decay widths (A15) and (A16) and read: and: T R,2 0.16γ 1 + 5 24γ 2 Notice that the decay rates (A15) and (A16), together with the inflaton decay width (30), also determine, via (A5), the fractional energy density of DR.

Inflationary observables and DM abundance
As in Sec. V, we start analyzing the cosmology of the model with two moduli by presenting the expression for the total number of e-foldings: where N reh is the duration of the reheating epoch after the end of inflation, while N φ1 and N φ2 denote respectively the number of e-foldings of the two EMD eras driven by the light moduli φ 1 and φ 2 , which look like: By rewriting N e in terms of fundamental parameters, we obtain: where we have set γ = 1. Notice that the total number of e-foldings does not depend on the initial misalignment value of the first modulus. As in the single modulus scenario, we obtain W 0 as a function of V by combining two constraints coming from the amplitude of the primordial scalar fluctuations and the geometrical requirement to have the volume of blowup modes hierarchically smaller than the overall internal volume. Following the same procedure as in Sec. V, we then extract the value of V from matching the observed DM abundance. Finally, this value of the volume fixes the DM mass for every combination of the underlying parameters. Interestingly, all data points correspond to a DM mass in the same range as in the single modulus case, m χ 10 10 -10 11 GeV, with a bias towards smaller values (65% of the data points result in m χ 10 10 GeV). In Fig. 7 we present the points in the (W 0 , V) plane which satisfy all our theoretical and phenomenological conditions. Note that 65% of the initial parameter set reproduces the correct DM abundance. Moreover we were not able to obtain values of n s within the 2σ range, as each point in Fig. 7 corresponds to n s at the lower end of the 3σ range. The qualitative behavior of the underlying parameters g s , Y φ2 , and λ is the same as in the single modulus case. Tab. II shows a representative choice of the parameters used to perform a numerical study of the full cosmological evolution in this scenario.   W0 and V, the resulting e-folding numbers, ns, mass scales, and inflaton coupling to hidden degrees of freedom c hid at a benchmark point that gives the right amplitude of the density perturbations and the correct DM abundance. The input parameters are gs = 0.1, Y φ 1 = Y φ 2 = 0.01, λ = 10 3 , N hid g = 12, and γ = 1.

Numerical analysis of cosmological evolution
We obtain the numerical evolution of energy densities for the scenario with two moduli. As in Sec. V B, we begin the evolution at H H inf with both σ and φ 1 oscillating. Though φ 2 begins oscillating shortly after this time, its energy density is subdominant and its initial evolution can be approximated as a matter component without altering the overall evolution. The other energy density components are again highly suppressed initially. The Boltzmann equations for this scenario are (neglecting the tiny production of DR from the decay of φ 1 ): dρ DR dt + 4Hρ DR = Γ σ→DR ρ σ + Γ φ2→DR ρ φ2 , dρ R dt + 4Hρ R = Γ σ→vis ρ σ + Γ φ1→vis ρ φ1 + Γ φ2→vis ρ φ2 , dn χ dt + 3Hn χ = Br χ Γ σ ρ σ m σ + σ ann v n 2 χ,eq − n 2 χ . (A26) A sample numerical solution is shown in Fig. 8 for the becnhmark point in Tab. II. The cosmological evolution is very similar to the one-modulus case. The effect of the second, lighter, modulus is to extend the EMD period to lower temperatures. Because the two moduli start with equal misalignments, we have a single extended period of EMD with a brief period of substantial radiation while the lighter modulus dominates, instead of two EMD periods separated by a RD phase. Fig. 9 shows the visible sector temperature as a function of the scale factor, where we have taken the temperature dependence of the relativistic degrees of freedom into account.