Mechanism for dark matter depopulation

Early decoupling of thermally produced dark matter particles due to feeble interactions with the surrounding plasma typically results in their excessive abundance. In this work we propose a simple mechanism for dark matter depopulation. It relies on a specific cosmological evolution under which dark matter particles become temporarily unstable and hence decay away reducing the overall abundance. The instability phase may be followed by an incomplete regeneration phase until the final abundance is established. We explicitly demonstrate this mechanism within a simple toy model of fermionic dark matter and discuss how it can be implemented in theoretically well motivated theories, such as the minimal supersymmetric standard model (MSSM) and for fermionic dark matter in the scotogenic model.


Introduction
Stable weakly interacting massive particles (WIMPs) are the most compelling dark matter (DM) candidates in several extensions of the standard model, notably in supersymmetric models. Perhaps the most attractive aspect of the WIMP paradigm is the fact that hypothetical ∼ 100 GeV stable neutral particles, interacting with the known standard model particles with the strength of weak interactions and being once in chemical equilibrium with them, would populate the universe in abundance comparable to the observed DM abundance. However, recent null results from direct detection experiments including LUX [1], PandaX-II [2] and XENON1T [3] put further constraints on the WIMP interaction cross section. This, along with the absence of any credible evidence for any new physics at the Large Hadron Collider (LHC), has cast doubt on the simplest WIMP hypothesis. Low WIMP cross section inevitably results in the earlier WIMP decoupling from the primordial thermal plasma and may lead to an excessive abundance of WIMP dark matter.
This situation is not atypical. Many well-motivated theoretical models predict over-abundant dark matter. For example, bino-like dark matter within the MSSM is over-abundant for most of the parameter area, except the cases when the bino is nearly degenerate in mass with the sleptons. Another example is fermionic dark matter in the scotogenic model [4], where neutrino masses and lepton-flavor-violating processes constrain DM interactions to be small and thus lead to the early decoupling of dark matter in excessive amount.
The above ramification of the early WIMP decoupling can be altered if the WIMP temporarily becomes unstable and its abundance is reduced through its decays. 1 This could happen if the symmetry that ensures the stability of the WIMP gets spontaneously broken at some early cosmological era and is restored at later times. 2 A similar idea has been recently put forward in Ref. [12]. The authors of this paper have not specified the production mechanism for dark matter and presented a rather contrived model which requires a specific mass arrangement between the different fields involved. On the contrary, we concentrate on thermal WIMP models. In addition, the instability phase may be followed by a regeneration phase, where DM is partially regenerated by the decay of fields which have been thermally produced during the instability phase. Our mechanism of depopulation can be employed in a variety of theoretical frameworks as is discussed below. Since we are primarily interested in a qualitative picture in this work, we rely on approximate analytic solutions as opposed to a comprehensive numerical study of the models.
The paper is organized as follows. In the next section we discuss the depopulation mechanism for a generic DM model due to the existence of a DM instability phase in the early universe. In section 3 we discuss the phase of DM regeneration and apply the proposed mechanism to a simple specific simple model of fermionic DM in section 4. In the last section 5 we briefly discuss how the depopulation mechanism can be implemented within the MSSM and the scotogenic model and draw our conclusions.

Dark matter depopulation: general consideration
Consider generic DM particles with mass m χ that undergo scatterings off the standard model particles in the primordial plasma at some high temperatures. The scattering cross section is assumed strong enough to keep DM particles in thermal and chemical equilibrium in the early universe. For the evolution of the cold DM abundance one typically assumes that DM particles become non-relativistic before their decoupling. If DM particles become non-relativistic after the decoupling DM is generically over-abundant. Since our scenario involves a phase where the DM abundance gets reduced the assumption of non-relativistic DM decoupling is not necessary. However, we must ensure that DM becomes non-relativistic during the depopulation phase to avoid DM repopulation due to inverse decays.
The DM abundance yield Y χ (x) = n χ /s (n χ and s are particle number and entropy densities, respectively), evolves according to the following Boltzmann equation: where x ≡ m χ /T , m χ being the DM particle mass, and In the above equation s χ = (2π 2 /45)g * m 3 χ and H χ = (π 2 g * /90) 1/2 m 2 χ M P are, respectively, the entropy density and the Hubble expansion rate both evaluated at T = m χ . M P ≈ 2.44 × 10 18 GeV is the reduced Planck mass and g * denotes the number of relativistic degrees of freedom in thermal equilibrium at a given temperature. For our purpose it is safe to ignore changes in g * with temperature. We also approximate the thermally-averaged cross section σv by its leading partial wave and thus regard λ as an x-independent constant. The power of x is given by n = 0 for s-wave scattering and n = 1 for p-wave scattering.
Initially, x O(1), the number density is approximately equal to the equilibrium number where g χ measures the degrees of freedom per DM particle, e.g. g χ = 7 4 for a Majorana fermion. Note, for relativistic particles Y (eq) χ stays constant and so is Y χ (x) (δ = 0) up until the DM particles become non-relativistic, x > 1. The evolution of non-relativistic DM density is more complicated. Namely, for non-relativistic DM particles Y χ (x) decreases exponentially together with the equilibrium yield as the universe cools down. At sufficiently low temperature x f (x f ≈ 20 for non-relativistic DM), the DM number density drops to the point when scatterings of DM particles become rare and DM is no longer able to sustain itself in the thermal bath. Close to this freeze-out temperature departure from the equilibrium yield becomes significant and Y χ ≈ δ Y (eq) χ (x) with x ∼ x f and: while for relativistic DM still δ = 0. Subsequently at x x f and x 1, the non-relativistic equilibrium yield becomes negligibly small and can be dropped from eq. (1). The analytic solution for the abundance yield at present time Y χ,0 ≈ Y χ (x → ∞) then reads as: 0.208(0.278) gχ g * for relativistic fermionic (bosonic) DM are defined by matching the solutions in two different regimes. Using the expression in eq. (4) for δ(x f ) and x f 1 for non-relativistic, eq. (5) we obtain an approximate solution for the DM yield and thus find Since typically λ 1, the DM abundance of relativistically decoupled DM is typically several orders of magnitude higher than that of non-relativistically decoupled DM. Now let us assume that there is a cosmological phase transition at x a to a phase where DM becomes unstable due to the spontaneous breaking of the symmetry which stabilizes DM, followed by the restoration of this symmetry at x b . In this instability phase x ∈ [x a , x b ] the DM particle is allowed to decay. As we have mentioned previously DM must be non-relativistic during at least part of this instability phase to ensure that inverse decays are suppressed and thus x b > 1. The Boltzmann equation then gets modified as: where n = 0 for s-wave annihilation, n = 1 for p-wave annihilation.
is the thermally-averaged DM particle decay width with the zero temperature decay width Γ χ and the modified Bessel functions K n (x). For non-relativistic DM it is approximately given by the zero temperature decay width, since As we are mostly interested in the region of parameter space when χ is (close-to) non-relativistic and to allow for an approximate analytic solution, we approximate it by the zero temperature decay width. If DM is relativistic during the initial stages of the instability phase, an approximate solution is obtained by neglecting DM decay before it becomes non-relativistic.
The last term describes the change in entropy S during the thermal evolution. The main contributions to entropy production are particle decoupling, but also a first order phase transition may generate additional entropy. If the DM yield is close to its equilibrium value, the change in entropy can be parameterized by an entropy dilution factor S(x)/S(x 0 ) which relates the DM yields as Y χ (x) = Y χ (x 0 )S(x 0 )/S(x). The same conclusion holds if the DM distribution is much larger than its equilibrium distribution, i.e. Y χ Y (eq) χ , and annihilations can be neglected compared to decays. In the most general case there is no analytic solution. In the following we assume that the instability phase occurs well above T = 1 GeV and thus we can neglect entropy production from particle decoupling, i.e. the number of relativistic degrees of freedom g * is approximately constant, and we approximate entropy production during a first order phase transition by an entropy dilution factor. Thus in the following discussion we neglect the last term in eq. (7), but remind the reader that for a first order phase transition entropy dilution factors have to be included in the final solution.
The effect of DM decays on its final abundance critically depends on when the instability phase took place. If it happened well before the standard freeze-out, i.e., x f x b , the subsequent scatterings can repopulate dark matter particles leading essentially to the standard abundance (6). However, a dramatic change in the final DM abundance is expected, when the instability phase follows the freeze-out (out of equilibrium decays), i.e. x f x a , or DM particles freeze-out during (x a < x f < x b ) or just after the instability phase (x f ∼ x a ). We separately consider the cases whether the equlibrium abundance of χ can be neglected or not.

Negligible equilibrium abundance
If the standard freeze-out precedes the instability phase, x f x a , and the decay rate is sufficiently small such that the equilibrium yield is always negligible, then one can ignore the equilibrium yield in eq. (7), which then becomes the Bernoulli equation and can be solved analytically (in terms of the incomplete Γ function): We have used an asymptotic expansion of the incomplete Γ-function in x to make the approximation in the last line. This approximation is accurate for all practical purposes. The 2-by-2 scattering cross section enters the solution (9) explicitly only through the matching condition, i.e. λ Y χ (x a ) ≈ (n + 1)x n+1 f , see eq. (6). We observe an exponential suppression factor which has clear intuitive explanation: the reduction of DM number density is essentially defined through the decay rate of DM particles per the universe expansion rate and the duration of the instability phase, ∆t ∝ ( . The present abundance yield can be obtained from eq. (5) by the substitution: in the absence of any regeneration phase as discussed in the next section.
Consider now the case where the standard freeze-out follows the instability phase, i.e. x f ∼ x b . In this case the solution to eq. (7) is given by This must be compared with the standard solution (4). Unless the decay rate Γ χ is extremely small, the first term in the denominator dominates and Y Hence, unlike scatterings, decays keep the DM abundance exponentially close to its equilibrium value even at the freeze-out temperature. The final yield, therefore, is obtained from eq. (5) by the substitution Y χ (x f ) → Y (eq) χ (x f ) and may be substantially smaller than the standard abundance in eq. (6). Furthermore, consider the case where freeze-out happens during the instability phase, i.e. x a < x f < x b . The solution for x a < x f is given by eq. (10), while for x f < x b the solution is given by (9). We combine these solutions by imposing the matching condition: The present day abundance then is again obtained from eq. (5) via the substitution:

Sizable equilibrium abundance
If inverse decays are fast enough to keep DM in equilibrium, we obtain an approximate analytic solution for the present day abundance from eq. (5) by replacing Finally, the intermediate regime where the inverse decays are relevant, but not fast enough to keep DM in equilibrium, requires a numerical solution.

Dark matter regeneration: general consideration
During the instability phase, x ∈ [x a , x b ], the scalar S which spontaneously breaks the symmetry has a vanishingly small thermal mass close to the critical temperature. Similarly there may be Goldstone bosons or massive gauge bosons from the breaking of additional continuous symmetries.
All these light degrees of freedom may be produced during the instability phase via decays of the DM candidate and scatterings with SM particles. The mass of these particles increases during the instability phase with decreasing temperature, because the VEV of the scalar S increases, S 2 ∝ 1 − x 2 a /x 2 , and thus depending on the second phase transition the light degrees of freedom may stay relativistic during the instability phase or become non-relativistic, if their thermal mass exceeds the temperature before the second phase transition.
Ultimately after the symmetry is restored at x b all light degrees of freedom of S are heavy and follow a Maxwell-Boltzmann distribution. Relativistic degrees of freedom of S which become non-relativistic at x b are far from their equilibrium distribution. We first focus on them and discuss degrees of freedom which become non-relativistic during the instability phase below. We assume that all relativistic degrees of freedom of S are in kinetic and chemical equilibrium following a Bose-Einstein distribution during the instability phase (it is sufficient if they are at the end). We assume that the scalars are still in kinetic equilibrium with the SM thermal bath after the phase transition and thus have the same temperature. This is satisfied if the relaxation rate Γ relax Γ coll /N coll (expressed in terms of the collision rate Γ coll and the number of collisions N coll ) exceeds the Hubble rate H Γ relax T 3m S Γ coll which we assume in the following. 3 Their number density does not change during the symmetry-restoring phase transition, hence we have: As . Particles which became non-relativistic during the instability phase follow a Maxwell-Boltzmann distribution both before and after the second phase transition. As their mass may change during the phase transition, but their number density remains the same during the phase transition, we find where m S denotes the mass of the scalar at the end of the instability phase and m S the zero temperature mass. Thus the scalar may also be over-abundant, ), directly after the phase transition, although its abundance is closer to the equilibrium abundance compared to above.
As the scalars are heavier than the DM mass and non-relativistic, they quickly annihilate to SM particles and decay to the DM candidate and a SM particle. In analogy to DM annihilation these processes are described by 4 where Γ S is the zero temperature scalar decay width which is a valid approximation in this case for the thermally averaged decay width, since S is non-relativistic. The annihilation cross section of S into SM particles is parameterized by 3 This condition is generally fulfilled, if the scalar interacts via electroweak interactions: The collision rate is Γ coll G 2 F T 5 for a (non-relativistic) particle S scattering with a relativistic particle in the SM thermal bath, and thus kinetic equilibrium is achieved for T 4 G −2 F mS/MP (40 MeV) 4 (mS/TeV). 4 Note that DM annihilations are neglected, because we assume that DM freeze-out occurs before the regeneration phase, i.e. x f < x b . Figure 1: A schematic illustration of the DM abundance where the standard freeze-out is followed by an instability phase where DM decays (T a > T > T b ) and regeneration of DM (T < T b ) until the final abundance Y ∞ χ is established.
Note that the DM may in general not be in chemical equilibrium, i.e. Y χ = Y (eq) χ . We again approximate the thermally-averaged cross-section by its leading partial wave. The power of x is given by m = 0 for s-wave scattering and m = 1 for p-wave scattering.
There are several relevant timescales: The phenomenology is different depending on the relative ordering of the timescales. As As inverse decays become important very quickly in the simple example model and thus a numerical solution is required, we do not give an explicit analytic solution, but only state that similar approximations as in the previous section may be used to obtain analytic solutions for some regions of parameter space. The DM depopulation scenario is schematically summarized on Fig. 1. This obviously only shows one scenario. DM may also freeze-out during the instability phase or while still being relativistic.

Simple fermionic model
We illustrate the mechanism described in the previous section within a simple toy model which contains an additional electroweak scalar doublet η and a Majorana fermion N apart from the SM particles. Both η and N are odd under an imposed Z 2 symmetry (Z 2 : η → −η and N → −N ), while SM particles are Z 2 -even. We largely follow the notation in Ref. [13] and refer the reader to this reference for more details. This model has direct applications to the scotogenic model [4] with fermionic DM and bino-like DM in the MSSM.
The most general Lagrangian involving extra fields is given by where we chose the Majorana mass term of the SM singlets to be diagonal with with real couplings λ i and µ 2 i and use the freedom of phase redefinitions of H and η to choose λ 5 < 0. We focus on the case with µ 2 η > 0 5 . We briefly summarize details on the different vacuum states, scalar masses, and stability conditions in App. A. It is convenient to define the parameter: Finite-temperature corrections to the potential can be parameterized by temperature-dependent masses to leading order in the high-temperature limit where we neglected all other Yukawa couplings apart from the top quark Yukawa coupling y t . From the stability conditions in eqs. (45) we can infer that c 1 + c 2 > 0. We also restrict our discussion to the parameter space without charge-breaking vacuum states and thus impose λ 4 < |λ 5 | and thus the VEVs are described by Then there are four charge-conserving minima (where we indicate the VEVs in brackets): (i) the electroweak conserving minimum EW c (v H = v η = 0); (ii) the phenomenologically desired zero temperature minimum I 1 (v H = 0, v η = 0), and two Z 2 -breaking minima, (iii)

Thermal evolution
The thermal evolution can be described by the trajectory in the plane (μ 2 H (T ),μ 2 η (T )), wherē We are interested in a thermal history which starts in the electroweak conserving phase EW c with µ 2 H (T ) < 0 andμ 2 η (T ) < 0, undergoes a phase transition to one of the phases I 2 /M such that µ 2 η (T ) > 0 and ends in I 1 withμ 2 H (T ) >μ 2 η (T ). During the intermediate phase the discrete Z 2 symmetry is broken by v η = 0 and finally restored after transitioning to I 1 . Following Ref. [13] we identify two distinct thermal evolutions of our interest that take place for different values of parameter R (19) which are shown in Fig. 2.
For R > 1 there are two phase transitions with an intermediate phase I 2 as shown in Fig. 2 (a). The first phase transition EW c → I 2 occurs at T EW c,2 while the second phase transition I 2 → I 1 happens at T 2,1 and thus the instability phase is given by the interval [x a , x b ] with: While the phase transition at x a is a second-order phase transition, the second phase transition at x a is a first-order phase transition (at x b the two minima are degenerate and separated by a barrier). We assume the phase transition to be instantaneous and estimate the entropy production ∆s = Q I 2 →I 1 /T b from the latent heat density Q I 2 →I 1 released during the phase transition [13] Q 5 The alternative scenario with µ 2 η < 0 requires large and negative H − η couplings which are in conflict with the stability conditions for the potential in this model.
The entropy dilution factor is then given by where g s * (T b ) denotes the number of relativistic degrees of freedom (as they enter the entropy density) just before the phase transition at temperature T b .
For 0 < R < 1 the phase diagram is more complicated as it contains three second-order phase transitions: EW c → I 2 → M → I 1 for 0 < R < 1, see Fig. 2 (b), with the corresponding critical temperatures T EW c,2 , T 2,M , and T M,1 , respectively. The instability phase is defined by:

Evolution of DM abundance
Having defined the above cosmological scenarios, we may proceed to calculate the DM abundance using the general formulae in sections 2 and 3. During the instability phase DM undergoes interactions through scattering and decay processes. The thermally averaged DM annihilation cross section is p-wave suppressed and can be approximated by the scalar masses can be neglected compared to the DM mass [14]. The DM decay (N → ηL a ) rate for heavy SM singlet fermions M N m W,Z is given by [15] Γ Note that there is an extra temperature dependence in the above equation which comes through thermal VEVs, v η , v H . This dependence can be neglected if v η v H and it drops out for the scenario with R > 1 since v H = 0 during the instability phase.
Next we turn to the regeneration phase, where η particles produced during the instability phase decay back to DM particles, η → N L a , with width Large M η − M N M La approximation was used in the above formula. The scalar mass M η weakly depends on the temperature T , but is well approximated by its zero temperature value and thus dominated by µ η . In addition, η annihilates into electroweak gauge bosons. In the limit of M η m W,Z and neglecting the mass splitting induced by the VEVs we obtain for s-wave annihilation cross section of the electroweak doublet scalar η the following decay width

Numerical Example
As a proof of principle we give one numerical example which leads to the correct DM relic density (Ω DM h 2 ) obs. = 0.112 ± 0.006 [16]. In particular, we examine the case explained in section 2.2 where the freeze out occurs during the phase transition such that the decays cause the dark matter to track the equilibrium yield until the phase transition ends. The quartic coupling λ 1 is fixed by the SM Higgs mass. In addition we choose the fermion mass M N , the charged scalar mass M η ± , Yukawa couplings and the quartic couplings λ 2,3,4,5 as follows All other parameters are chosen at their respective best-fit values [17].

Discussion and conclusions
In this paper we have proposed a mechanism of DM depopulation, which occurs due to the temporary breaking of a symmetry that ensures DM stability. During the instability phase the DM particles decay and their relative abundance drops down. Through these decays the primordial plasma gets also populated by additional species, which subsequently (after the instability phase) decay back to DM particles. The final DM abundance is established after this incomplete regeneration.
We have demonstrated the DM depopulation mechanism using a simple extension of the standard model which contains a DM fermion N and an additional electroweak doublet scalar field η, both of which are Z 2 -odd. The latter field is assumed to be an inert field, i.e., without VEV at zero temperature. However, it plays a crucial role cosmologically in setting up the instability phase with broken Z 2 symmetry. With a suitable choice of parameters one can reduce an initially-large DM abundance down to the observed value.
The fermionic DM model we have explicitly discussed has no deep theoretical motivation, but rather serves an illustrative purpose. The mechanism of depopulation can be applied to various theoretically better motivated models. In fact, our model is the closest analogue to the scotogenic model [4] which incorporates neutrino masses together with fermionic DM within the SM extension with two electroweak doublet scalars and three SM singlet fermions. Fermionic DM which is thermally produced via annihilations is in conflict with search for lepton-flavor-violating processes [18]. Conversely, satisfying those constraints leads to low annihilation rates and thus overproduction of DM particles. The mechanism of depopulation can be straightforwardly applied to this class of models, resolving the tension with observations.
Another class of well-motivated models with multiple scalar states are supersymmetric theories. The neutral lightest supersymmetric particle (LSP) (typically the neutralino) may serve as a natural candidate for cold DM in models with R-parity conservation. However, for significantly large space of parameters in supersymmetric extensions of the SM (especially for bino-like neutralinos) theory predicts DM to be over-abundant. The mechanism of depopulation can be applied to such models as well. Let us outline here one of the possible implementations of our mechanism within the minimal supersymmetric SM, where R-parity is broken temporarily during the instability phase of the sneutrino condensate. Imagine the relevant light states at low energies are just the SM Higgs boson and one sneutrino (the rest of sparticles are assumed to be heavy).
The zero temperature scalar potential in this limit reads: where mν is the sneutrino soft supersymmetry breaking mass parameter. Since we are interested in a phase where the sneutrino develops a non-zero vacuum expectation value, we allow it to be tachyonic 6 . The physical sneutrino mass, however, receives contribution from the electroweak symmetry breaking and radiative corrections. The mass square, M 2 ν = m 2 ν + m Z m h 2 + rad. corr. , must obviously be positive, hence the following constraint must be met: The dominant high temperature corrections to the potential (34) reads (we ignore linear and logarithmic in T terms): The analysis of the full potential V 0 + V T reveals that at sufficiently large T the Higgs and the sneutrino fields reside in the origin, h T = ν T = 0, and hence electroweak symmetry and R-parity are exact. Once the universe cools down to Tν c ≈ 2.78|mν|, an instability develops inν-direction and the sneutrino field develops a non-zero vacuum expectation value, ν T = v m Z −m 2 ν − ανT 2 1/2 , while the Higgs field remains at the origin. Further cooling down to T h c ≈ 143 GeV results in the electroweak phase transition due to the Higgs field condensate h T = v m h m 2 h − 2α h T 1/2 . Positive contribution to the sneutrino mass parameter from the Higgs condensate starts to dominate and brings the sneutrino field back to the origin, restoring R-parity. Hence, for suitable m 2 ν we indeed account for a phase in the early universe, where R-parity is broken spontaneously. During this phase the neutralino LSP ceases to be a stable particle. More specifically, condensation of the sneutrino field leads to a spontaneous breaking of R-parity and to a mixing of neutralinos and neutrinos. Through this mixing the neutralino LSP decays into standard model particles, the dominant decay channel being 2-body χ → Zν for m χ > m Z +m ν . The longitudinal degrees of freedom of Z-boson during the instability state become massive sneutrino states and decay back to neutralino DM during the regeneration phase. Although the full phenomenological validity of this particular supersymmetric model requires further study, we are confident that the mechanism of depopulation can be implemented within supersymmetric models along the lines outlined here.
In conclusion, we find that the proposed depopulation mechanism can be implemented within various DM models to bring the DM abundance to observed valued, which otherwise would be considered empirically invalid.
We do not report the masses for the vacuum M where both VEVs, v H and v η , are non-zero, because the expressions for the masses in the vacuum M are more complicated and we not refer to them in the main part of the manuscript. After electroweak symmetry breaking the masses of the electroweak vector bosons are given by

A.2 Relevant cross sections and decay widths
The DM annihilation cross section of fermionic DM N into leptons is p-wave suppressed. Thus to leading order the annihilation cross section is given by [14] vσ(N N → νν, + − ) = v 2 r 2 (1 − 2r + 2r 2 ) 24πM 2 with the relative velocity v. The thermally-averaged DM annihilation cross section is For heavy DM M N M η R , M η I during the instability phase, the DM annihilation cross section can be approximated by vσ(N N → νν, The decay rate of DM to leptons is [15] Γ N = Γ W + − + Γ W − + + Γ Zν (55) with the leptonic mixing matrix U and the active-sterile mixing ξ N Y vη √ 2M N for vanishing final state lepton masses. Note that the Fermi constant in the Z 2 -breaking phase is determined by . Summing over all flavors in the final state we find for M N m W,Z