Neutralino dark matter in scenarios with early matter domination

We investigate the production of neutralino dark matter in a cosmological scenario featuring an early matter dominated era ending at a relatively low reheating temperature. In such scenarios different production mechanisms of weakly interacting massive particles (WIMPs), besides the well-studied thermal production, can be important. This opens up new regions of parameter space where the lightest neutralino, as the best-known supersymmetric (SUSY) WIMP, obtains the required relic abundance. Many of these new sets of parameters are also compatible with current limits from colliders as well as direct and indirect WIMP searches. In particular, in standard cosmology bino-like neutralinos, which emerge naturally as lightest neutralino in many models, can have the desired relic density only in some finetuned regions of parameter space where the effective annihilation cross section is enhanced by co-annihilation or an s-channel pole. In contrast, if the energy density of the universe was dominated by long-lived PeV-scale particles (e.g. moduli or Polonyi fields), bino-like neutralinos can obtain the required relic density over wide regions of supersymmetric parameter space. We identify the interesting ranges of mass and decay properties of the heavy long-lived particles, carefully treating the evolution of the temperature of the thermal background.


Introduction
The lightest neutralino as lightest supersymmetric particle (LSP) is one of the oldest and most studied examples of a weakly interacting massive particle (WIMP) candidate for the cosmological Dark Matter (DM); see e.g. [1] for an early exploration of parameter space, and [2,3] for reviews. The minimal supersymmetric extension of the Standard Model (MSSM) contains four neutralino current eigenstates: a bino, a wino, and two higgsinos. Given current collider constraints on superparticles, in particular on the masses of charginos and the heavier neutralinos, we now know that over most of parameter space, the mass eigenstates are relatively pure states, with little mixing.
Most analyses of WIMP DM worked in the framework of standard cosmology, where the Universe was radiation-dominated starting at the end of inflation and ending at a temperature around 1 eV. Moreover, it is usually assumed that the post-inflationary reheat temperature was sufficiently high that WIMPs attained full thermal (chemical and kinetic) equilibrium. The WIMP relic density is then basically inversely proportional to its (effective) annihilation cross section [4,5]. In that case higgsino-like WIMPs typically need to have a mass near 1 TeV to have the correct relic density, and a wino-like WIMP should be at least two times heavier. While it has recently been pointed out that these values might be lowered by 30% or so due to co-annihilation effects [6], the required values are still uncomfortably high when compared to estimates of weak-scale finetuning in the MSSM. In particular, while bounds on the masses of scalar tops and gluinos based on simple loop JHEP12(2018)042 calculations [7] are somewhat controversial [8][9][10], it is generally agreed that higgsino, and hence LSP, masses above several hundred GeV would lead to percent level (or worse) finetuning; note that in the MSSM the higgsino mass enters the relevant finetuning condition already at tree-level. 1 In standard cosmology, higgsino-or wino-like WIMPs with masses in the few hundred GeV range would have too small a relic density. In contrast, a bino-like WIMP has too large a relic density in such a scenario, unless its effective annihilation cross section is boosted by co-annihilation [12][13][14] or by an s-channel pole [12,15].
A predicted underdensity of WIMP DM can be cured by adding another DM component, e.g. axions [16][17][18]; this can be done within the framework of minimal cosmology, and without changing TeV-scale particle physics. On the other hand, a scenario that predicts too large a relic density for a given DM candidate is clearly excluded. This argument thus disfavors bino-like WIMPs, at least within minimal cosmology.
At the same time bino-like WIMPs quite easily satisfy the increasingly stringent constraints from direct WIMP searches [19,20]; these searches exclude many scenarios where the WIMP is higgsino-like, if the latter contributes most or all of DM. Moreover, indirect searches [19,21] now exclude models where most or all of DM consists of wino-like (higgsino-like) WIMPs with mass below ∼ 0.8 (∼ 0.4) TeV, but hardly constrain the parameter space if the LSP is bino-like. These null results therefore favor bino-like WIMPs. At the same time bino-like neutralinos often emerge as LSP in simple models where the superparticle spectrum can be described by a small number of free parameters. In particular, if gaugino masses unify at or near the same scale where the gauge couplings meet in the MSSM, the weak-scale bino mass will be about half of the wino mass. Moreover, if stop squarks and Higgs bosons have similar soft breaking masses at this very high energy scale, the weak-scale higgsino mass parameter typically comes out larger than the bino mass.
These arguments motivate us to investigate a non-minimal cosmological scenario, in the hope of finding an extended region of parameter space where a bino-like WIMP obtains the required relic density. In particular, we analyse scenarios featuring an early matterdominated epoch sometime between the end of inflation and Big Bang nucleosynthesis (BBN). This is quite well motivated, since UV-complete theories like supergravity [22] and superstring theory often contain heavy but long-lived scalar particles, nowadays usually called moduli. They are long-lived since their couplings to MSSM fields are suppressed by the inverse of the Planck mass. Nevertheless they can attain large densities if their mass is below the Hubble parameter during inflation [23][24][25][26][27]. The success of standard BBN implies that the moduli-dominated epoch should end at a final reheat temperature of at least 4 MeV [28][29][30][31][32].
It has been pointed out more than ten years ago that the WIMP relic density in this scenario can be either smaller or larger than in standard cosmology [33,34]. On the one hand, since the moduli decay out of equilibrium, they increase the entropy density of the universe, thereby diluting a pre-existing WIMP density. At the same time new WIMP production mechanisms become possible, including direct moduli to WIMP decays. In

JHEP12(2018)042
addition to the effective WIMP annihilation cross section and mass, the final WIMP relic density now also depends on the mass and lifetime of the modulus as well as on the effective branching ratio for modulus to WIMP decays.
The impact of an early matter dominated epoch on WIMP DM has been studied before [33][34][35][36][37][38][39][40][41][42][43]; specifically supersymmetric WIMPs were considered in this context in [44][45][46][47][48]. However, as we pointed out recently [49], an accurate treatment of the radiation component during the decay of the moduli is very important: even though the entropy is no longer conserved in this scenario, one still has to normalize the WIMP number density to the known radiation density, which is related to the entropy density ; any inaccuracies in the calculation of the latter therefore immediately affect the final prediction of the WIMP relic density. Note also that the temperature of the thermal background during modulus domination can vary over several orders of magnitude. An accurate treatment of the temperature dependence of the effective number of degrees of freedom therefore becomes mandatory if one aims for precise predictions. Here we use the treatment of ref. [50], which in turn is based on lattice QCD analyses of the equation of state, which determines the relation between temperature and entropy and energy density in QCD. Moreover, we use the latest results from direct and indirect WIMP searches to further constrain the allowed parameter space. We find that large regions of parameter space with good bino-like WIMP indeed survive, if moduli are not too heavy and have a sufficiently small branching ratio for decays into WIMPs.
The remainder of this article is organized as follows. In section 2 we discuss the set of equations we need to solve in order to compute the WIMP relic density. Then section 3 briefly reviews neutralino DM in standard cosmology, including a discussion of current experimental constraints. In section 4 we return to moduli cosmology, delineating regions in parameter space where various WIMP production mechanisms are important. We then present numerical results from a random scan over the MSSM parameter space. Finally, we summarize our results in section 5.

Basic framework
As stated in the Introduction, we use the formalism of ref. [49], which carefully treats the evolution of the radiation component with temperature, using [50] for the precise evolution of the effective number of degrees of freedom with temperature.
In this article we ignore the details of thermalization process [51][52][53] by working in the approximation that the decay products of the massive "modulus" particle φ (other than the WIMP χ) thermalize instantaneously. In particular, we do not consider the production of WIMPs from reactions involving energetic φ decay products prior to their thermalization, which can be important for WIMP masses not too much smaller than the φ mass [54]. In this approximation we only need to consider the integrated Boltzmann equations (for densities in normal space, rather than for phase space densities) describing the evolution of the mass density ρ φ of the massive, very weakly interacting particle φ, of the entropy JHEP12(2018)042 density s R of the radiation component, and of the WIMP number density n χ : The first eq. (2.1) describes the reduction of the φ density through Hubble expansion and φ decay, Γ φ being the total φ decay width which is the inverse of the lifetime τ φ . The second equation describes the evolution of the thermal background, whose entropy density is diluted by the Hubble expansion but increased by φ decay into radiation as well as by out-of-equilibrium annihilation of WIMPs; however, this last term is always negligible in practice. The entropy production from φ decay depends on the effective branching ratiō χ the average number of χ particles produced per φ decay and M φ the mass of φ. Since we are interested in scenarios with M χ M φ ,B 1 in all cases of interest. 2 Finally, the last eq. (2.1) describes the evolution of the WIMP number density, which is again diluted by Hubble expansion and enhanced by φ → χ decays; the last term of this equation describes annihilation of WIMPs into SM particles as well as inverse annihilation of SM particles, assumed to be in thermal equilibrium, into WIMPs.
The Hubble parameter appearing in eqs. (2.1) is as usual determined by the Friedmann equation: where M Pl 1.22 · 10 19 GeV is the Planck mass. Since the χ number density becomes comoving constant well within the radiation-dominated epoch after φ decay, the last term in eq. (2.3) is very small for the range of temperatures over which we need to solve the Boltzmann equations numerically, but the contribution from radiation is important; it is given by

4)
g * being the number of relativistic degrees of freedom defined via ρ R . The temperature T is computed from the entropy density, which is given by

JHEP12(2018)042
where h * is another measure for the number of relativistic degrees of freedom. A massless, non-interacting particle contributes equally to g * and h * , but particle masses as well as the strong interactions around the QCD deconfinement transition affect g * and h * differently [50]. As mentioned in the Introduction, the relic density of the lightest neutralino is often affected by co-annihilation effects. In particular, if the LSP is wino-like, its mass is very similar to that of the lightest chargino; a higgsino-like LSP is close in mass to both the lightest chargino and the second-lightest neutralino. We include co-annihilation using the framework developed in ref. [12]. Note that reactions turning the LSP into one of its co-annihilation partners or vice versa remain in equilibrium well after the system of superparticles decouples from the SM particles. Moreover, after this decoupling eventually each superparticle will decay into an LSP [plus some irrelevant SM particle(s)]. These observations allow to interpret n χ appearing in eqs. (2.1) as the sum over the number densities of all (relevant) superparticles, n χ = i n i , and n eq = i n eq,i for the equilibrium densities. The effective thermally averaged cross section for superparticle annihilation can then be written as [46,55] σv eff = i=1 j=1 σ ij v ij n eq,i n eq n eq,j n eq . (2.6) Here the double sum runs over all superparticles that can co-annihilate with the LSP. In practice it is sufficient to include weakly interacting particles whose masses are within about 20% of the LSP mass, whereas strongly interacting superparticles within 30 or 35% of the LSP mass should be considered. The modulus dominated epoch ends when most φ particles have decayed. Here we are interested in scenarios with a fairly long early matter dominated era. Today's radiation then almost all comes from φ decays. In the instantaneous decay approximation, the "reheat" temperature after φ decay is given by (2.7) In practice this equation has to be solved iteratively, due to the temperature dependence of g * appearing on the right-hand side (rhs). Since we treat φ decays exactly, rather than in the instantaneous decay approximation, T RH strictly speaking has no physical meaning; it nevertheless remains a good measure characterizing the end of the early matter dominated epoch. The only free parameter appearing in eq. (2.7) is the φ decay width. We make the usual assumption that its coupling to (some) (MS)SM particles are suppressed by a single power of the Planck mass, in which case The value of the constant C is model dependent. In the rest of this paper we assume α = 1. In scenarios with sufficiently long early matter domination, the final result depends essentially only on Γ φ and B χ /M φ once all densities have been normalized to today's JHEP12(2018)042 radiation density. Results for α = 1 can thus be read off from our numerical results by choosing a different value of M φ such that Γ φ is the same as in our ansatz, and then rescaling B χ so that the original value of B χ /M φ is restored; this will work as long as the resulting M φ M χ and the resulting B χ < 1. For our numerical work we follow ref. [40] and introduce dimensionless quantities: where a is the scale factor in the Friedmann-Robertson-Walker metric. The Hubble parameter can then be written as where c 1 = 3 8π and the dimensionless Hubble parameterH is given by We also use the dimensionless quantity A to parameterize the evolution of the Universe, rather than the (cosmological) time t appearing in the original Boltzmann equations (2.1). The new, dimensionless evolution equations are thus: 3 (2.12) Here we have introduced c ρ = π 2 g * (T RH ) 30 . We solve these differential equations with the initial conditions [37,41,49] The case with non-vanishing initial WIMP and radiation densities has been studied in [49]. There we found that an initial R I ≤ Φ I does not change the result if γ 10 20 . The following argument shows that this is not difficult to achieve. The primordial modulus field, created during inflation, begins to behave like an ensemble of free particles at H = M φ . The earliest plausible H I is thus given by M φ . H I 10 20 Γ φ therefore requires M φ 10 20 Γ φ , which via eq. (2.8) implies M φ 10 −10 M Pl . On the other hand, the requirement that

JHEP12(2018)042
WIMPs do not attain thermal equilibrium after φ decay requires T RH M χ /20, and hence via eq. (2.7), M φ 8 · 10 7 GeV [M χ /(1 TeV)] 2/3 , which is a significantly stronger bound for M χ ≤ 1 TeV. Moreover, at H = M φ the universe will only have been radiationdominated if at that time the temperature exceeded 2 · 10 12 GeV M φ /(10 7 GeV) 1/2 . If the temperature at H = M φ was below this value, the universe was φ matter dominated for all H between M φ and Γ φ , and the initial radiation content can be safely ignored; and even if the temperature at H = M φ was somewhat above this (already quite high) value, the modulus dominated epoch might still have been long enough to make the initial radiation, and χ, content irrelevant. As mentioned above, the final DM relic density is obtained by normalizing to the known radiation density [49]: (2.14) Here Ω i denotes today's mass or energy density of component i in units of the critical density (i.e., Ω tot = 1 yields a flat universe), h is the Hubble constant in units of 100 km/(Mpc · s), and T now is the current temperature of the cosmic microwave background (CMB) radiation. A E is the value of A where the numerical solution of the Boltzmann equations is terminated. It should be well within the radiation-dominated epoch after all φ particles have decayed and χ has completely decoupled, so that R and X become constant. Note that the entropy density is comoving-constant for A > A E . Finally, we use the PDG values for today's quantities [19]: The present DM relic density is also quite well known [19] 4 Since we normalize to today's radiation density, the final result is essentially independent of the initial φ density Φ I , as long as the period of early matter domination lasted sufficiently long (γ 10 10 ) and R I , X I 1. The reason is that in this limit all densities are proportional to Φ I , which therefore drops out in the ratio used in eq. (2.14).

JHEP12(2018)042
In the MSSM with R-parity conservation [56,57] the lightest supersymmetric particle is stable; it can be a neutralino which is therefore a potential DM candidate [1]. The MSSM contains four neutralino interaction (or current) eigenstates: a bino, a wino, and two higgsinos, where the former are superpartners of the neutral U(1) Y and SU(2) gauge bosons and the latter are superpartners of the two MSSM Higgs doublets. Upon electroweak symmetry breaking these states mix via a 4 × 4 mass matrix. However, over most of parameter space the mixing is rather small, i.e. usually the four mass eigenstates are dominated by a single current eigenstate, with small admixtures of the three other states. This is true because searches at LEP and the LHC require [19] the wino and higgsino mass parameters to be significantly above the mass of the W bosons, which sets the scale for the off-diagonal entries in the neutralino mass matrix. The neutralino masses are therefore essentially set by the soft breaking bino mass M 1 , the soft breaking wino mass M 2 and the supersymmetric higgsino mass µ. Note that the bino is a gauge singlet. M 1 is therefore still essentially unconstrained experimentally, if one does not assume it to be related to M 2 [58].
If |M 1 | < |M 2 | , |µ| the LSP will be bino-like. In most of parameter space its annihilation cross section is dominated by the exchange of sfermions in the t− and u-channel [15], so that where the last factor is due to the P -wave suppression of the annihilation of a Majorana fermion into massless SM fermions; here α em is the fine structure constant and θ W is the weak mixing angle. Recall that in standard cosmology, the WIMP relic density is inversely proportional to its (effective) annihilation cross section [4]. The increasing lower bounds on the sfermion masses [19] imply that the cross section (3.1) is too small, leading to bino overdensity in standard cosmology. This can only be evaded if the bino co-annihilates, e.g. with aτ slepton [13] or at squark [14], or if the annihilation cross section is enhanced by a nearby resonance, in particular for M 1 m A /2 where m A is the mass of the neutral CP-odd Higgs boson of the MSSM [15] For |µ| < |M 1 | , |M 2 | the LSP is higgsino-like and can annihilate efficiently into W + W − and Z 0 Z 0 pairs, via the exchange of other higgsino-like chargino and neutralino states, respectively. In this case annihilation from an S-wave initial state is allowed, leading to a thermally averaged cross section [59] σH v = g 4 512πµ 2 21 + 3 tan 2 θ W + 11 tan 4 θ W . (3. 2) The presence of these other higgsino-like states, whose masses are within a few GeV of that of the LSP, implies that co-annihilation effects are always important for higgsino-like LSP [55,60]. Here also other final states contribute (e.g. W ± Z 0 forχ ± 1χ 0 1 co-annihilation, or ff pairs via Z 0 exchange forχ 0 1χ 0 2 co-annihilation, where f is an SM fermion). All these cross sections also scale like 1/µ 2 , so that in standard cosmology Note that Sommerfeld enhancement is quite small for higgsino-like LSP [61].

JHEP12(2018)042
Finally, for |M 2 | < |M 1 | , |µ| the LSP is wino-like, with thermally averaged cross section [59] σW v = 3g 4 16πM 2 2 . (3.4) As in the case of higgsino-like LSP, co-annihilation is always important for wino-like LSP, since the charged wino is very close in mass to the neutral one. The effective cross section determining the wino relic density in standard cosmology therefore differs somewhat from eq. (3.4), but it still satisfies Even in the absence of Sommerfeld enhancement [62,63] a wino-mass near 2.5 TeV is required in order to obtain the desired relic density (2.16) in standard cosmology. Note also that a wino-like LSP is not compatible with the unification of gaugino masses near the scale where the MSSM gauge couplings unify; in such models M 2 2M 1 at the weak scale.
In standard cosmology there are no moduli fields, and the comoving entropy density is conserved during WIMP decoupling. We therefore only have to solve the second eq. (2.12). We used MicrOMEGAs [64] for this purpose; 5 this program also allows to compute the LSP annihilation cross section in today's universe, needed to check indirect detection constraints, and the LSP-nucleon scattering cross sections, which are constrained by direct detection experiments. Moreover, we used SuSpect [65] for the computation of the spectrum of superparticles and Higgs bosons in the MSSM. The scans over parameter space where performed with the help of the T3PS [66] which parallelizes the computation, leading to great gains in speed. This is not crucial for standard cosmology, but becomes important in the non-standard scenario, where the solution of the Boltzmann equations becomes considerably more computationally intensive.
We performed two different scans. The first is in the framework of the phenomenological MSSM (pMSSM), where the soft breaking parameters are defined directly at the TeV scale. Here we follow ref. [46] and parameterize the spectrum with 10 independent free parameters ("p10MSSM"): the masses of the three MSSM gauginos M 1 , M 2 , M 3 ; the trilinear scalar soft breaking parameters fort squarks andτ sleptons, A t , A τ ; the higgsino mass µ; the mass m A of the CP-odd neutral Higgs boson; a common soft breaking mass m Q 3 for third generation squarks; a common soft breaking mass m L 3 for third generation sleptons; and the ratio of vacuum expectation values tan β. Since the A parameters are multiplied with the corresponding Yukawa couplings, they are irrelevant for the first and second generation, so we set them to zero. The A parameter forb squarks has also been fixed; note that L − R mixing in theb sector is dominated by µ, not by A b , since µ comes with a factor tan β here. The sfermion masses of the first and second generation are chosen to lie sufficiently far above M 1 that these sfermions do not contribute significantly to co-annihilation; however, co-annihilation with third generation sfermions is possible. 5 In MicrOMEGAs the thermally averaged cross section is a function of temperature T . Even in nonstandard cosmology WIMPs will be in thermal equilibrium until the temperature falls well below the WIMP mass. This function is thus only needed for T Mχ. In our numerical calculation we therefore set σv (T ≥ Mχ/13) = σv (T = Mχ/13), keeping the full temperature dependence only for T < Mχ/13. This helps to speed up the calculation without significant loss of accuracy in the final result.
We vary all these parameters simultaneously in a random scan, using flat distributions for the values between the limits shown in table 1. We only consider points where the predicted mass of the lighter CP-even Higgs boson lies within 3 GeV of the measured value of 125 GeV [67,68]. Our scan probably includes some points with first and second generation squark or gluino masses below current bounds [19]; however, the precise values of these parameters basically do not affect the physics of neutralino DM, apart from possible co-annihilation (which our lower limit of the squark mass range excludes, as mentioned above). The lower limit of the scan range of third generation squark masses is also quite low; however, the requirement of a sufficiently heavy Higgs boson requires TeV-scale stop masses, which satisfy current LHC search limits. Since in the pMSSM scan M 1 , M 2 and µ are varied independently, all three types of LSP can occur.
We also have performed another scan over the parameter space of the constrained MSSM (cMSSM), where the soft breaking parameters are assumed to unify near the scale of Grand Unification, M X 2·10 16 GeV. Specifically, at this very high scale all gauginos are assumed to have the same mass m 1/2 , and all scalars (sfermions as well as Higgs bosons) get a common soft breaking mass m 0 . The trilinear soft breaking terms also unify to A 0 . tan β is again a free parameter. We chose the sign of µ to be positive, as in the pMSSM scan; this has little effect on the relic density, but slightly increases the spin-independent scattering cross section on nucleons for bino-like LSP [69]. The physical masses of superparticles and Higgs bosons are again computed with the help of SuSpect, which also solves the relevant renormalization group equations. The range of parameters we used is shown in table 2.
In the cMSSM the LSP is usually bino-like. Since gaugino masses are assumed to unify, a wino-like LSP is not possible in this scenario. Moreover, the large top Yukawa

JHEP12(2018)042
Parameter Range scalar mass 0.1 < m 0 < 6 gaugino mass 0.1 < m 1/2 < 6 trilinear coupling −12 < A 0 < 12 ratio of Higgs doublet VEVs 1 < tan β < 60 sign of µ parameter µ > 0 Table 2. The range of cMSSM parameters we used in our scan are given at the GUT scale where the MSSM gauge couplings meet. All dimensionful parameters are given in TeV. Examples of the predicted thermal relic density are shown in figure 1, for the p10MSSM (left) and cMSSM (right). In the left frame we observe two bands of points with relic density ∝ M 2 χ . Here the LSP is higgsino-like (green points) or wino-like (blue points); the observed behavior conforms with the expectations of eqs. (3.3) and (3.5). In contrast, the red points, where the LSP is bino-like, are widely scattered, since the effective bino annihilation cross section depends not only on the LSP mass, but also on sfermion masses; in addition, schannel Higgs exchange resonances can play a role. As already noted in the Introduction, scenarios with bino-like LSP typically lead to (much) too large a relic density, whereas higgsino-or wino-like LSPs attain the correct relic density only for masses that require quite severe electroweak finetuning.
As expected, there are no points with wino-like LSP in the cMSSM scan (right frame of figure 1), and very few points where the LSP is higgsino-like. Moreover, the number JHEP12(2018)042 of points with rather light bino-like LSP is also relatively small. This is because large gaugino, and hence bino, masses help to increase the weak-scale stop masses; as noted above, quite large stop masses are needed in order to reproduce the Higgs mass of about 125 GeV. Furthermore the relic density of bino-like LSPs tends to be even higher than in the pMSSM scan, due to the reduced opportunities for co-annihilation and resonant s-channel enhancement, as explained above.
We now turn to the direct detection constraints depicted in figure 2. Note that this constraint, as well as the constraint from indirect detection discussed below, assume that the given LSP forms all of dark matter (in our galaxy). We saw that this is frequently not the case in standard cosmology. Since these constraints do not directly depend on the cosmological scenario, we nevertheless already discuss them here.
We see that most of the (red) points with bino-like LSP are safely below the bound. In the left, pMSSM, frame the few red points above the bound have quite small masses for first generation squark and/or the heavy neutral Higgs boson; such scenarios are very difficult to realize in the cMSSM, so in that scenario nearly all points with bino-like LSP survive this constraint. In contrast, many points with higgsino-like LSP are excluded even in the pMSSM case. In the cMSSM scan all points with higgsino-like LSP are excluded. The reason is that in the cMSSM even for the points with higgsino-like LSP the higgsino mass is not much smaller than that of the bino. This leads to relatively large higgsinobino mixing, which yields relatively large couplings of the LSP to neutral Higgs bosons, and hence rather large scattering cross sections. We finally note that most points with wino-like LSP in the pMSSM scan survive this constraint even if winos form all of DM.
Finally, the thermally averaged neutralino annihilation cross section in today's universe is depicted in figure 3. Since WIMPs are now very non-relativistic, this is essentially the S-wave contribution in the limit v → 0. The solid lines are exclusion contours from a combination of recent γ telescope data, for different assumptions of the dominant final JHEP12(2018)042 state. These bounds have again been obtained under the assumption that a given WIMP forms all of DM.
We see that these constraints exclude wino-like neutralinos with M χ ≤ 0.8 TeV, and higgsino-like LSP with M χ ≤ 0.4 TeV. Note that co-annihilation does not play any role in today's universe. Wino-and higgsino-like LSPs both annihilate predominantly into heavy gauge bosons, i.e. the relevant exclusion limit is the one shown by the cyan curves. In contrast, all scenarios with bino-like LSP are allowed. The safety margin is even higher in the cMSSM scan, since a large s-channel resonance enhancement of the annihilation cross section is even less likely there, as noted above. Figures 2 and 3 clearly favor bino-like neutralinos, especially for the "natural" region of parameter space where M χ 0.5 TeV; higgsino-like states near that mass are also acceptable, if somewhat marginal in view of finetuning and approaching constraints from both direct and indirect searches. On the other hand, figure 1 shows that in standard cosmology, such binos usually have much too high a relic density, whereas the relic density of 0.5 TeV higgsinos is too low. This motivates the investigation of our non-standard cosmological scenario, to which we now return.

Neutralino production in non-standard cosmology
As noted above, an early matter dominated epoch can either increase or reduce the predicted WIMP relic density, depending on the values of the free parameters. The most relevant ones are the ratio of the "reheating" temperature T RH at the end of that epoch, given by eq. (2.7); the mass M χ of the WIMP; the WIMP annihilation cross section, which (by crossing) also determines the rate of WIMP production (or "inverse annihilation"); and the combination M χ B χ /M φ , which determines the importance of direct φ → χ decays. In this section we first make some general remarks on the parameter space, and then present

JHEP12(2018)042
some examples of scenarios leading to the approximately correct relic density for bino-or higgsino-like neutralinos with "natural" masses.

Discussion of the parameter space
Here we analyse the dependence of various LSP production mechanisms on the mass M φ , which determines the reheat temperature via eq. (2.7), and on the branching ratio for φ → χ decays.

Modulus mass
Let us first discuss the relevant range of M φ . Its lower limit is set by the requirement that modulus decay should not interfere with big bang nucleosynthesis. This requires [32] T RH ≥ 4 MeV, and hence via eqs. (2.7) and (2.8), On the other hand, if φ is very heavy, T RH becomes so large that WIMPs will maintain (or reach) full thermal equilibrium even after all φ particles have decayed. This would bring us back to standard cosmology as far as WIMP physics is concerned. Recalling that in standard cosmology, WIMPs decouple at T = T F M χ /20, and taking into account that in reality the end of φ domination is not sharp, we define the upper end of the interesting range of M φ through the condition T RH ≤ M χ /10, which implies In both eqs. (4.1) and (4.2) we have again assumed that the effective φ coupling parameter α = 1. Eq. (2.8) implies that both bounds on M φ should be multiplied by α −1/3 if α = 1.
It is important to note that even for φ masses well below the upper bound (4.2) thermal contributions to the final LSP relic density can be important. This is because φ decays quickly create a thermal background [37], whose maximal temperature exceeds T RH by a factor γ 1/4 [49]. If the effective LSP annihilation cross section is not too small, LSPs will therefore typically attain full thermal equilibrium during the φ matter dominated epoch. The resulting contribution to the final LSP relic density scales like [37] where x FO = Mχ T FO , T FO being the LSP decoupling temperature in the φ dominated epoch. Note that this contribution to the relic density again scales like the inverse of the effective LSP annihilation cross section, but the coefficient is smaller than in standard cosmology by a factor ∝ (T RH /T FO ) 3 . For a given LSP mass and annihilation cross section, T FO is somewhat larger than the freeze-out temperature in standard cosmology, since for a given temperature the Hubble parameter is larger by a factor (T /T RH ) 2 ; this translates in a stronger M χ dependence of x FO ≡ M χ /T FO . Moreover, a decrease of T RH implies an increase of T FO , reducing the suppression factor (T RH /T FO ) 3 even further.

JHEP12(2018)042
For small M φ , and hence small T RH , and very small annihilation cross section, LSPs may not attain thermal equilibrium at all during the φ matter dominated epoch. In case of neutralino LSPs, this may happen for M φ 10 6 GeV and M χ 1 TeV. Even in this case there will in general be a thermal contribution to the final LSP relic density, due to the production of superparticles from the thermal plasma. As long as the LSP annihilation cross section is below that required to obtain full thermal equilibrium, this contribution will be proportional to this cross section [37]. Finally, independent of whether LSPs attain thermal equilibrium during φ domination or not, the thermal contribution to the LSP relic density is bounded from above [37]: (4.4) For our choice α = 1 this means that thermal LSP production can become relevant only for M φ 3 · 10 6 GeV (100 GeV/M χ ) 8/15 . Note that the upper bound (4.4) can get saturated only for rather small annihilation cross section; for higgsino-or wino-like neutralinos, which attained full thermal equilibrium during φ domination for all cases of interest, the maximal thermal contribution is significantly smaller.

Branching fraction
In addition to thermal production, neutralinos can also be produced from φ decays. Note that here not only direct decays into the LSP are relevant, but decays into all superparticles, since the heavier superparticles will quickly decay into the lightest neutralino. Moreover, most φ particles decay at a temperature near T RH . The upper bound (4.2) ensures that χ particles are not in thermal equilibrium at that temperature. 6 The final LSP relic density can then be written as [41] Ω χ h 2 = Ω χ, thermal h 2 + Ω χ, decay h 2 , (4.5) where Ω χ, thermal h 2 originates from χ interactions with the thermal plasma and has been discussed above. Here we analyze the second contribution, which is due to φ → χ decays. Although by assumption χ particles were not in thermal equilibrium at T T RH , the χ density may attain "quasi-static" equilibrium, where the r.h.s. of the second eq. (2.12) vanishes due to a cancellation between the term ∝ X 2 and the term ∝ B χ . This will happen if the LSP annihilation cross section is above a critical value [41], which scales like M φ /(B χ M Pl T 2 RH ). Using eqs. (2.7) and (2.8) and writing σv = κ/M 2 χ where κ is dimensionless, this can be translated into a critical value of B χ : In the second step we have given approximate numerical values for higgsino-and wino-like LSP. The annihilation cross section of bino-like LSPs does not simply scale like M −2 χ times JHEP12(2018)042 a constant; in fact, in most cases of interest the bino annihilation cross section is below the critical value. If the annihilation cross section is well above its critical value or, equivalently, if B χ > B χ, crit , Ω χ h 2 is independent of B χ and scales inversely with the effective annihilation cross section [41]. The dependence on σv is thus the same as in standard cosmology, but the relic density is parametrically enhanced by a factor T F /T RH , where T F M χ /20 is the LSP decoupling temperature in standard cosmology. Schematically, (4.7) In the opposite limit, where the annihilation cross section is well below its critical value, Ω χ, decay h 2 is simply proportional to B χ : The factor T RH /M φ results because the φ number density n φ = ρ φ /M φ , and ρ φ is related to the entropy density s after φ decay via ρ φ ∝ T RH s, with s being removed by the final normalization to today's photon density. In the second step we have used T RH ∝ M 3/2 φ , see eqs. (2.7) and (2.8).

Numerical results
We are now ready to present some numerical results yielding cosmologically interesting relic densities for neutralinos, in particular for relatively light bino-or higgsino-like neutralinos; we saw in section 3 that indirect searches already exclude wino-like neutralinos as main contributor to the total DM density for "natural" masses, which we interpret as M χ 500 GeV. We use the same MSSM parameters as in section 3. The resulting relic density can therefore be directly compared to figure 1. We use the same computer codes as for the case of standard cosmology, including some custom-written new routines for MicrOMEGAs that solve the coupled system of Boltzmann equations (2.12).
We structure this discussion by analyzing different values for M φ within the range defined by eqs. (4.1) and (4.2).

Light moduli
We first discuss the case M φ = 5 · 10 5 GeV, corresponding to T RH = 40.68 MeV; recall that we assume α = 1 in eq. (2.8). This scenario is rather simple to analyse. On the one hand, the bound (4.4) implies that neutralino production from the thermal bath is negligible in this case. 7 Moreover, eq. (4.6) shows that the annihilation of neutralinos produced in φ decay become significant only for B χ 10 −3 even for wino-like LSP. These considerations are confirmed by figure 4. In the left frame B χ = 10 −5 . In that case neutralino annihilation is always negligible, hence the result is independent of σv . On the other hand, for B χ = 10 −3 (right frame) the annihilation of wino-like neutralinos becomes relevant for M χ < 1 TeV, but it is still basically irrelevant for bino-and higgsinolike neutralinos.

JHEP12(2018)042
This allows to determine the first choice of parameters giving the desired relic density for "natural" bino-or higgsino-like LSP: This equation is quite accurate for M φ = 5 · 10 5 GeV for which the results of figure 4 have been obtained. It also works fairly well for other values in the range 10 5 ≤ M φ ≤ 10 6 GeV, where the lower bound had been given in eq. (4.1). For M φ 10 6 GeV, T RH comes close to or exceeds the QCD transition temperature, so that the effective numbers of degree of freedom g * and h * show a sizable dependence on T RH and hence on M φ , which is not captured by the simple expression (4.9). For slightly larger M φ contributions to the final LSP relic density from interactions with the thermal plasma can also become significant.
Nevertheless eq. (4.9) remains a reasonable approximation even for M φ 10 7 GeV as long as the annihilation cross section is sufficiently small, e.g. for bino-or higgsino-like neutralinos with M χ 300 GeV. Moreover, it works for all versions of the MSSM, and in fact for all WIMPs whose annihilation cross section does not substantially exceed that of wino-like neutralinos, simply because WIMP annihilation is negligible in these cases. The only difference in the cMSSM is that a wino-like LSP cannot be realized, and a higgsino-like LSP is rare, as already discussed in section 3. We therefore do not show numerical results for the cMSSM in this light moduli scenario.

Intermediate-mass moduli
We now turn to the case M φ = 5 · 10 6 GeV, corresponding to T RH = 848.60 MeV. This case is more complicated than the previous scenario, since now thermal and non-thermal LSP production can both be important; annihilation of non-thermally produced neutralinos also plays a prominent role for some range of parameters leading to a cosmologically interesting DM relic density. This is illustrated in figure 5. The top-left frame is for B χ = 10 −3 . We see that this relatively large branching ratio leads to a sizable overdensity if the neutralino is bino-like. The fact that most of the red points lie along a line implies that bino annihilation is still insignificant in most cases. This implies that the branching ratio (4.9) can be interpreted as an upper bound for intermediate values of M φ , if the LSP is bino-like. 8 On the other hand, the annihilation cross sections for higgsino-and wino-like neutralinos with M χ 1 TeV are already so large that eq. (4.7) applies, where the final DM

JHEP12(2018)042
relic density is independent of B χ as long as B χ > B χ, crit . At the same time, the thermal production of higgsino-and wino-like LSPs still leads to a very small contribution to the final LSP relic density. Note that for M χ ∼ 170 GeV the green points in this frame lie near the desired relic density. This allows to identify a second region of parameter space that will give the correct DM density, this time for higgsino-like LSPs: The bound on B χ is a numerical approximation of the requirement B χ > B χ, crit , see eq. (4.6), for the range of higgsino masses leading to approximately the correct DM relic density. Of course, M φ should not be so large that higgsinos thermalize at T T RH , i.e. the bound (4.2) should also hold for eq. (4.10) to be applicable.
Comparing this figure with the left frame of figure 1 we see that the gap between the blue and green bands is even somewhat larger here than in standard cosmology. As discussed earlier in this section, in both cases the relic density scales like 1/ σv . However, in standard cosmology it also scales linearly with x F = M χ /T F , which grows logarithmically with increasing annihilation cross section. This factor, which does not exist in eq. (4.7), slightly reduces the difference of the predicted relic density of wino-and higgsino-like neutralinos in standard cosmology.
In the top-right frame we have reduced B χ to 10 −5 . This lies below B χ, crit for M χ 1 TeV in all cases. In fact, eq. (4.9) is still applicable in this region of parameter space, but the required LSP mass near 1 TeV is already outside the range we consider to be natural.
For slightly smaller M χ the green and blue dots begin to diverge. Here LSP annihilation after φ decay begins to be relevant, although even at M χ 100 GeV the ratio between the green and blue points is smaller than in the top-left frame, i.e. we have not fully reached the regime described by eq. (4.7).
Moreover, bino-like LSPs, which have a much smaller annihilation cross section, now begin to receive sizable contributions from the thermal plasma. For this value of M φ neutralinos with M χ 1 TeV will always attain thermal equilibrium rather early in the φ matter dominated epoch. For wino-or higgsino-like neutralinos this contribution to the final relic density is still much smaller than the non-thermal contribution from φ decays, even after late neutralino annihilation is included. However, for bino-like neutralinos with very small annihilation cross sections this thermal contribution can reproduce the required relic density! This leads to another region of parameter space with the required relic density, this time for bino-like states: This is a rough approximation. For example, we have ignored the factor x 4 FO in eq. (4.3); recall from the discussion of that equation that x FO has a stronger dependence on M χ than the corresponding quantity in standard cosmology does. Moreover, we have expressed the thermally averaged cross section in GeV units, with 10 −13 GeV −2 being near the smallest cross section we found for bino-like LSP. An increase of the annihilation cross section can JHEP12(2018)042 Figure 6. The neutralino relic density vs neutralino mass for M φ = 5 · 10 7 GeV and B χ = 10 −3 (10 −9 ) in the left (right) frame. The notation is as in figure 4. be compensated by a slight decrease of MB, or by an even smaller (relative) increase of M φ . We should emphasize that eq. (4.11) is applicable only if B χ is well below the value given in (4.9). One can also find combinations of parameters where φ →B decays and thermalB production both give comparable contributions to the final relic density; these can be obtained quite easily from eqs. (4.9) and (4.11).
In the bottom-left frame of figure 5 we have reduced B χ even further, to 10 −7 . This is well below B χ, crit for all neutralinos, but the resulting non-thermal contribution is still much larger than the thermal contribution for higgsino-or wino-like neutralinos. On the other hand, the final relic density of bino-like neutralinos is now dominated by the thermal contribution for M χ 400 GeV. This contribution drops quite quickly with increasing M χ . This is mostly from the explicit M −3 χ factor in eq. (4.3); the fact that the cross section also tends to increase with M χ for bino-like LSP, as shown in figure 1, also contributes. Finally, x FO decreases logarithmically with increasing M χ . As expected, the red points near Ω χ h 2 = 0.1 that appeared in the top-right frame also show up here, since for these points the non-thermal contribution was negligible already in the former case, and is even smaller here.
The bottom-right frame of figure 5 is again for B χ = 10 −7 , but this time shows results for the cMSSM scan. The red points here show a similar trend as in the bottom-left frame, showing that the region approximated by eq. (4.11) can be accessed in the cMSSM as well.

Heavy moduli
Finally, we consider a scenario with M φ = 5 · 10 7 GeV, near the upper end (4.2) of the range where a φ dominated epoch can modify the final LSP relic density. In fact, since now T RH = 25.49 GeV, we expect that the neutralino relic density should essentially coincide with that in standard cosmology if M χ 300 GeV. This is borne out by figure 6. In the left frame we chose B χ = 10 −3 . This is well above B χ, crit for both higgsino-and wino-like neutralinos. Their relic density scales like the inverse of the cross section, and hence essentially like M 2 χ . The slight change of slope of the JHEP12(2018)042 lines of green and blue points around M χ 300 GeV occurs because in standard cosmology, which is applicable for M χ below this point, the explicit x F factor in the expression for the relic density slightly softens the dependence on the cross section, and hence on M χ , as remarked in the discussion of the top-left frame of figure 5. Note that thanks to this steepening of the slope, higgsino-like neutralinos obtain the required relic density already at M χ 0.6 TeV here, well below, and hence more natural than, the value near 1 TeV required in standard cosmology. In fact, in this region of parameter space eq. (4.10) is still approximately applicable.
In contrast, nearly all red points indicating a bino-like neutralino lie well above the required value, the exception being scenarios where the effective annihilation cross section is enhanced by co-annihilations or s-channel resonances. As expected, for M χ 300 GeV the distribution of red points is very similar to that in the left frame of figure 1, which use the same choices of MSSM parameters. The accumulation of red points with M χ 700 GeV results from φ →B decays, where lateB annihilation is still negligible, even for these large densities, which are more than a factor of 100 above the desired value.
In the right frame of figure 6 we therefore took a very small branching ratio, B χ = 10 −9 . The non-thermal contribution to the neutralino relic density is then always negligible. As expected, for M χ 300 GeV the green and blue bands are the same as in the left frame. However, for M χ 1 TeV the relic density is essentially set by thermal freeze-out during the φ matter dominated epoch. Eq. (4.3) shows that the relic density there scales like 1/(M 3 χ σv ) ∝ 1/M χ , which explains the fall-off with increasing M χ at these large neutralino masses. In between there is therefore a region where the relic density depends only very weakly on the LSP mass, even though LSP production is purely thermal here. Unfortunately the value of the relic density in this region of parameter space is about one order of magnitude below the required value even for higgsino-like neutralino.
Most of the parameter points with bino-like LSP still predict a significantly too large relic density. However, for M χ 300 GeV the upper envelope of the region populated by the red points decreases quickly with increasing M χ ; this is the same behavior we saw in the two lower frames of figure 5. As a result, for M χ 500 GeV the density of red points with Ω χ h 2 0.12 is significantly higher than in standard cosmology, cf. figure 1. The upper envelope even approaches the desired value for M χ 2.5 TeV. However, since these points violate our naturalness criterion, we do not attempt to define another interesting region of parameter space in which these points lie.

Summary and conclusions
In this paper we investigated supersymmetric neutralino dark matter in the framework of a non-standard cosmological scenario with an early matter dominated epoch. Building on our earlier work [49], which improved the accuracy of the solution of the relevant Boltzmann equations through a careful treatment of the thermal medium, we looked for regions of parameter space where relatively light neutralinos can form all of dark matter; our focus on neutralinos with mass at or below 500 GeV is motivated by naturalness arguments.

JHEP12(2018)042
After setting up the basic framework, in section 3 we reviewed neutralino DM within standard cosmology. In agreement with many earlier studies, we found that a bino-like neutralino typically has too high a relic density, whereas higgsino-or wino-like neutralinos can obtain the desired relic density only for masses well above the "natural" range. Moreover, under the assumption that neutralinos form all of DM, indirect searches exclude wino-like DM for masses below about 0.8 TeV, which is already well above our naturalness cut-off. For higgsino-like neutralinos the corresponding bound is around 400 GeV, leaving some range of masses where higgsino DM could be (barely) natural if it could get the required relic density. In contrast, if the lightest neutralino is bino-like neither direct nor indirect DM searches are very constraining. 9 In section 4 we therefore set out to find regions of parameter space where a light bino, or a higgsino with mass near 400 GeV, can obtain the required relic density. The main new parameters, in addition to those already present in standard cosmology, are the mass M φ of the heavy particle φ that accounts for the early matter-dominated epoch, and its branching ratio B χ into the DM candidate. We found three distinct regions, described by eqs. (4.9), (4.10) and (4.11); this is the main result of the present paper. In particular, for relatively small M φ , below 10 6 GeV, neutralino annihilation is negligible, and a simple relation for the required B χ as a function of M φ and M χ results, see eq. (4.9); this works for both bino-and higgsino-like states. For higgsino-like neutralinos, annihilation becomes important for M φ > 10 6 GeV, leading to a second range where φ →H decays followed byH annihilation produces the desired relic density; this is still a purely non-thermal mechanism. Finally, for M φ 5 · 10 6 GeV a third region opens up, where bino-like neutralinos obtain the desired relic density due to thermal freeze-out during the early matter dominated epoch, see eq. (4.11).
In this analysis we treated M φ and B χ as completely free parameters, and assumed that the φ decay width scales like M 3 φ /M 2 Pl . In particular, in agreement with more general results found in ref. [49], we found that bino-like neutralinos can only have the desired relic density even in this non-minimal cosmological scenario if B χ 10 −4 . It would be important to find concrete models possessing a φ particle with the desired properties. In particular, the upper bound on B χ may not be easy to satisfy once higher-order decays into three-or even four-body final states have been included. We leave this investigation to future work.
T kd is well above TRH. This is because the period of early matter domination will lead to an enhanced growth of structure at very small length scales [72]. These "microhalos" will be destroyed by free streaming of DM particles if T kd TRH. In case of WIMPs, T kd TRH > 4 MeV requires sfermion masses of tens of TeV, well above the range probed in our scan.

JHEP12(2018)042
Forschungsgemeinschaft via the TR33 "The Dark Universe". FH is supported financially by the Deutsche Akademische Austauschdienst (DAAD). FH also thanks the hospitality and support of Goethe University of Frankfurt at the final stage of this work.
Note added. While completing this work, ref. [74] appeared, which has some overlap with our work. Their results for a cosmology with early matter domination qualitatively agree with ours, but the main focus of their analysis is different.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.