Dark Matter and Baryogenesis from Visible-Sector Long-Lived Particles

We present a minimal extension of the standard model that includes a long-lived fermion with weak-scale mass and an ${\cal O}({\rm GeV})$ fermionic dark matter candidate both of which are coupled to quarks. Decays of a TeV-scale colored scalar in a radiation-dominated phase bring the former to a thermal abundance while also producing dark matter. The long-lived fermion then dominates the energy density of the Universe and drives a period of early matter domination. It decays to reheat the Universe, mainly through baryon-number-violating interactions that also generate a baryon asymmetry, with a small branching fraction to dark matter. We find the allowed parameter space of the model and show that it can be probed by proposed long-lived particle searches as well as next-generation neutron-antineutron oscillation experiments. This model provides a robust explanation of dark matter and baryogenesis as long as the Universe is in a radiation-dominated phase at $T \gtrsim {\cal O}({\rm TeV})$.


I. INTRODUCTION
There are various lines of evidence pointing to the existence of dark matter (DM) in the Universe [1]. However, the identity of DM remains a profound problem at the interface of cosmology and particle physics. In addition, explaining the observed DM abundance depends on the details of the thermal history of the early Universe. Thermal freeze-out in a radiation-dominated (RD) Universe can yield the correct relic abundance for a specific value of the (thermally averaged) annihilation rate σ ann v = 3 × 10 −26 cm 3 s −1 . On the other hand, in nonstandard thermal histories this can be achieved for both much larger or smaller values of σ ann v [2]. Important classes of early-Universe models (notably those arising from string theory) typically lead to nonstandard thermal histories that involve one or more epochs of early matter domination (EMD) [3]. An EMD phase is driven by a matter-like species that comes to dominate the energy density of the Universe at early times and subsequently decays to establish RD prior to big bang nucleosynthesis (BBN). This can be due to coherent oscillations of a scalar field like a string modulus displaced from the minimum of its potential during inflation, or from long-lived nonrelativistic quanta produced in the postinflationary Universe [4][5][6][7]. Scenarios with EMD have novel predictions that can be tested in observations and indirect DM searches [8][9][10][11][12].
We recently proposed a scenario [13] where a visible-sector long-lived particle (LLP) with a weak-scale mass drives an EMD epoch. The LLP reaches thermal equilibrium at temperatures well above its mass, maintains a frozen comoving number density, and dominates the energy density of the Universe. It eventually decays to standard-model (SM) particles which can occur all the way up to the onset of BBN. It was shown that this scenario can accommodate the observed DM abundance for both large and small values of σ ann v . Moreover, the parameter space of this scenario may be probed by the proposed LLP searches at the Large Hadron Collider (LHC) [14][15][16].
Here, we present an explicit realization of this scenario that yields the correct DM abundance and gives rise to successful baryogenesis. The model is a minimal extension of the SM that involves three new fields with baryonnumber-violating couplings to quarks: an O(TeV) colored scalar X and two Majorana fermions N (a neutral LLP) and χ (the DM candidate) with masses m N ∼ 100 GeV and m χ ≈ O(GeV) respectively. N acquires a thermalequilibrium abundance through X decays at T m N and, due to its long lifetime, dominates the energy density of the Universe and drives an epoch of EMD. It eventually decays to establish a RD Universe before the onset of BBN. The main contributions to the DM relic abundance come from X and N decays, while baryon-number-violating decays of N also generate a baryon asymmetry through one-loop electroweak corrections.
We find the allowed parameter space of the model where the observed DM abundance and correct baryon asymmetry are obtained for some benchmark points that are within or close to the LHC reach. The allowed region overlaps with the sweet spot of the proposed searches for hadronically decaying LLPs [14][15][16]. Moreover, the model also has predictions for baryon-number-violating processes at low energies. As we show, its parameter space can be probed by next-generation neutron-antineutron oscillation experiments. A positive signal in these experiments, combined with the requirement to satisfy very stringent bounds from double proton decay, will then help us further constrain the baryon-number-violating couplings of the LLP in this model. This model is robust as it is largely independent of the details of the postinflationary history provided that the Universe is in a RD phase at T O(TeV).
The rest of this paper is organized as follows. In Section II, we present the model and obtain the resulting DM relic abundance and baryon asymmetry. In Section III, we present our main results and identify the allowed region of the parameter space of the model. In Section IV, we discuss experimental signals of the model within this allowed region focusing on LLP searches and low-energy processes. We conclude the paper in Section V. Some details of our calculations are presented in the Appendix.

II. THE MODEL
The model is based on a minimal extension of the SM that was proposed for low-scale baryogenesis and DM [17] (see [18] for a supersymmetric version). Using two-component Weyl fermions, the Lagrangian is: where u c and d c denote the left-handed up-type and down-type antiquarks respectively, and h ij is an antisymmetric tensor. Flavor indices are denoted by i and j, while color indices are omitted for simplicity. X is an iso-singlet color-triplet scalar of hypercharge +4/3. N and χ are singlet fermions 1 with m χ m N m X . As discussed in [17], χ is absolutely stable if m p − m e ≤ m χ ≤ m p + m e . It is therefore a DM candidate in this mass window and, in addition, can help address the baryon-DM coincidence puzzle.
We denote the maximum value of |h i |, |h ij |, |h i | by h, h , h respectively, and quarks and antiquarks by q andq. Then, the (rest frame) decay width of X is given by: where the factor of 2 in front of h 2 accounts for the two color combinations inqq. The rate for N self-annihilation at energies E m X follows: where n N is the number density of N . The factor of 3 accounts for the color of q in the final state. N can also annihilate with quarks, which results in: where n q is the quark number density, and the factor of 6 is due to the possible color combinations in qqq. The rest-frame width for three-body decay of N is: Note that this decay violates baryon number and, as we will describe below, can be the origin of baryon asymmetry in the model. There is also two-body radiative decay of N at the one-loop level. The corresponding decay width, after all multiplicity factors for the particles in the loop are taken into account, is given by [19]: where α em is the electromagnetic fine structure constant. The total decay width of N is This model leads to the following thermal history, starting with a RD phase at an initial temperature T i m X , with H being the Hubble rate: (1) H H(T = m X ). Initially, X is in thermal equilibrium due to its gauge interactions with the SM particles. Its decay and inverse decay to N and χ then results in (for details, see Appendix A): where F is given in Eq. (31) and n eq N,χ denotes the equilibrium value of the N or χ number density. At H ∼ H(T = m N ), N quanta become nonrelativistic and their comoving energy density becomes frozen, remaining constant for H H dom provided that: which ensures that N self-annihilation and its annihilation with SM particles are inefficient at T m N . N starts to dominate the energy density of the Universe once its energy density becomes comparable to that of radiation. This happens at H = H dom , which is found to be [13]: where g * dom denotes the number of relativistic degrees of freedom at N dominance, and g * N is taken at T = m N 3 . A necessary condition for N dominance is: otherwise, N will decay before dominating and the Universe will remain in a RD phase.
(3) H dom > H Γ N . The Universe enters an EMD era that lasts until N quanta decay, at which point the Universe is reheated to a temperature T dec given by T dec = (90/g * dec π 2 ) 1/4 (Γ N M P ) 1/2 . This must happen before the onset of BBN 4 , which translates to T dec 7 − 8 MeV in our case. The entropy released by N decay dilutes any preexisting relic abundance by a factor d, which is given by [13]: where T dom = (45/g * dom π 2 ) 1/4 (H dom M P ) 1/2 . We see from Eqs. (11) and (13) that: The normalized number density of N at the time of its decay, called the "yield" Y N , is given by: We see that F (γ N ) drops out of the final expression for Y N , which can be understood as a smaller n N implies a later N dominance and hence a smaller dilution factor 5 . The relic density of DM in our model is given by: where: 2 More precisely, g * n N ∝ a −3 as long as N is in chemical equilibrium with the thermal bath, while n N ∝ a −3 when it is chemically decoupled. The difference is negligible though as the number of relativistic degrees of freedom, g * , changes minimally for m N T < m X . 3 To compute the number of relativistic degrees of freedom at a given temperature, g * (T ), we have made use of the smooth function generated from the data presented in Table S2 of [20] with cubic spline interpolation. 4 An O(1) MeV lower bound on the final reheat temperature is studied in [21] considering the thermalization process for neutrinos including neutrino oscillations and self interactions for both radiative and hadronic decays. A bound of the same order is obtained in [22] from CMB anisotropies considering detailed production of the relic neutrino background with three-flavor oscillations. 5 We note that if the condition in Eq. (12) is not satisfied, there will be no EMD epoch. In this case, we have d = 1 and Y N < 3T dec /4m N .
The first term on the right-hand (RH) side of Eq. (16) is the contribution from X decay 6 , while the second term accounts for production from N decay. This must match the observed relic abundance that follows: where the last relation is due to m χ ≈ m p . Given that F (γ N ) ≤ 1, we see from Eqs. (14) and (16) that: Matching the observed DM abundance thus requires that T dec 10 −9 m N if F (γ χ ) ∼ 1 or Br N →χ ∼ 1. The lower limit of T dec 7 − 8 MeV from BBN then results in m N 10 7 GeV (and hence m X 10 7 GeV). This is far beyond the reach of the LHC. Therefore, to have any hope of testing the model in the proposed LLP searches, we must have F (γ χ ) 1 and Br N →χ 1. These imply that DM production from the thermal bath occurs in the freeze-in regime, and N decay to final states with DM is subdominant.
The baryon-number-violating interactions of X in Eq. (1) allow for this model to generate the observed baryon asymmetry of the Universe. Baryogenesis from two-body decays of X (if m X > m N ) or N (if m N > m X ) has been studied in the non-supersymmteric [17,23] and sypersymmetric versions of the model [18,[24][25][26][27]. The latter scenario cannot work in our case as we have taken N to be lighter than X. Also, out-of-equilibrium decay of X will yield too little asymmetry because of Boltzmann suppression of n X at T < m X and subsequent dilution by the EMD phase.
However, as shown in [18], three-body decays of N can also be the source of baryogenesis. The baryon asymmetry of the Universe thus produced is given by: The asymmetry parameter B is determined by one-loop electroweak corrections giving rise to the following typical leading term [18]: Here α 2 is the SU (2) W gauge fine structure constant, m W is the mass of the W boson, and subscripts in the numerator denote quarks.

III. RESULTS
In this section we present our results. We choose two benchmark points m X = 3 TeV with m N = 100 GeV, and m X = 1 TeV with m N = 100 GeV. This overlaps with the most important physics target of the MATHUSLA (MAssive Timing Hodoscope for Ultra Stable neutraL pArticles) detector concept [14], namely hadronically decaying LLPs with masses in the O(10 GeV) − O(100 GeV) range [16] (more on this later). It also puts X within or close to the LHC reach.
In Fig. 1, we show the region (shaded green) that corresponds to EMD driven by N in the h − h plane for our benchmark points m X = 3 TeV, m N = 100 GeV (left panel) and m X = 1 TeV, m N = 100 GeV (right panel). Below the boundary of the shaded region given by the solid yellow line, decay of N occurs after the onset of BBN. Above the solid blue boundary line, N decays before it has had a chance to dominate the energy density of the Universe. Furthermore, above the solid light-blue line, N self annihilation is efficient such that its abundance is significantly depleted. We additionally show a solid red line above which annihilations of N with the SM bath become efficient, thus the abundance of N is not depleted in the white band between the solid blue and red lines. Above the dashed black line, we have Γ X→N H(T = m X ), implying that F (γ N ) = 1 in Eq. (11) throughout the bulk of the green region. Below this line, F (γ N ) < 1 and N does not reach an equilibrium abundance.
The dotted black curves depict contours of h that achieve the observed DM relic abundance. These contours have two main behaviors depending on the source of the dominant contribution to the DM abundance: h ∝ 1/h when thermal production from X decay in the bath is dominant, while h ∝ h when direct production from decay of N is dominant 7 . We see that the thermal contribution is dominant in the majority of the parameter space shown in the figure. Note that DM is produced in the freeze-in regime as we have F (γ χ ) 1 as long as h 10 −7 for m X ∼ TeV. Additionally, for the contribution from N decay, two-body decay dominates over three-body decay by approximately a factor of 0.15 (m X /m N ) 2 , independent of the couplings as can be seen from Eqs. (5) and (6).
We see that successful baryogenesis can be accomplished in the shaded parameter space for our benchmark masses. Due to the smallness of h , N dominantly decays to three-body final states qqq andqqq throughout the allowed parameter space. For the chosen value of m N = 100 GeV, Eq. (21) gives B ∼ 3 × 10 −8 and, in fact, the maximum value of B is obtained for m N ∼ 100 GeV 8 . The dashed magenta line in Fig. 1 corresponds to the observed value of the baryon asymmetry η B ∼ 10 −10 and it sits within the green EMD region close to the blue boundary which corresponds to moderate dilution. For the small values of h in the figure, this line scales with m X in the same way as the blue boundary line, however B introduces additional dependence on m N . Thus for higher values of m X and m N the baryogenesis line would be pushed outside of the shaded region toward the top-right. In the region between the red and blue lines, where the N abundance is not depleted by annihilations but there is no period of EMD, the observed DM abundance can be obtained for h ≈ 6 × 10 −11 . However, successful baryogenesis would require h ≈ 25h which would therefore overproduce DM in this region.
The effect of very energetic DM particles on the matter power spectrum can further constrain the parameter space. If DM particles are highly relativistic at the time of kinetic decoupling from the thermal bath, their free streaming can prevent the formation of structure at small scales such as those observed in the Lyman-α forest and Milky Way satellite galaxies. As shown in Appendix C, this effect does not lead to additional constraints on the parameter space for our benchmark points, but can become important for heavier masses.

IV. EXPERIMENTAL SIGNALS OF THE MODEL
In this section, we discuss the prospects of probing this model via particle physics experiments. LLP searches -The long lifetime of N provides us with an opportunity to directly probe it via the proposed LLP searches at the LHC. The decay length l N of N particles produced from X decay in colliders is related to the rest-frame lifetime τ N via l N =bcτ N , whereb is the average boost factor of N . We find that l N > 100 m within the allowed parameter space, which makes N a natural target for dedicated LLP searches. In particular, MATHUSLA would search for neutral LLPs produced in the high luminosity LHC (HL-LHC) collisions, and could discover LLPs with decay lengths up to 3 × 10 7 m (corresponding to lifetimes close to the age of the Universe at the onset of BBN) [15]. We note that throughout the allowed parameter space of Fig. 1, three-body final states qqq andqqq totally dominate N decay. The lifetime of N in this region varies between 10 −6 s (blue boundary) and 3 × 10 −2 s (yellow boundary) 9 , and the corresponding decay lengths are: Low-energy processes -It is also possible to probe this model via low-energy experiments. The most stringent experimental bound on the model parameters comes from double proton decay pp → K + K + . This process is due to the ∆B = 2 and ∆s = 2 effective interaction (uds) 2 . The N -mediated and χ-mediated contributions to this operator are given by (h 1 h 12 ) 2 /16π 2 m 4 X m N and (h 1 h 12 ) 2 /16π 2 m 4 X m χ respectively [28], where we have included flavor indices. The current experimental limit of τ pp→K + K + > 1.7 × 10 32 yr [29,30], is then translated into an upper bound on h 1 h 12 for given values of m X and m N . Scaling the limit in [28] for our benchmark points results in: Another important limit is from n −n oscillation, which can be used to constrain h 1 h 13 2 and h 1 h 13 2 . The ∆B = 2 effective operator (udb) 2 gives rise to n −n oscillation at the one-loop level [28]. The contributions from N and χ couplings to the oscillation amplitude G n−n can be parameterized as h 2 1 h 13 4 m N /16π 2 m 6 X and h 2 1 h 13 4 m χ /16π 2 m 6 X respectively (times logarithmic factors ln(m 2 X /m 2 N ) and ln(m 2 X /m 2 χ ) respectively) [23]. The oscillation lifetime is given by τ n−n ∼ (Λ 6 QCD G n−n ) −1 . The current experimental limit is τ n−n ≥ 3 × 10 8 s [31][32][33], while the next generation experiments predict a sensitivity τ n−n ≥ 5 × 10 10 s [34]. For our benchmark points, and choosing Λ QCD = 250 MeV, we find the following range that is compatible with current limits and also within the reach of future experiments: We note that the values of h (which is the largest of h ij ) in Fig. 1 are so small that the upper limit on h 1 h 12 in Eq. (23) is comfortably satisfied and h 1 h 13 2 is well outside the range in Eq. (24). In Fig. 2 we show the flavordependent bounds from Eqs. (23,24) on hh and hh 2 . In the left-panel, corresponding to m X = 3 TeV and m N = 100 GeV, the double proton bound (violet dash-dot line) is inside the green region, and the n −n oscillation band (shaded orange) from Eq. (24) overlaps with a portion of the green region. In the right panel, where m X = 1 TeV and m N = 100 GeV, the situation is similar with slightly more overlap between the orange and green bands.
Given that h and h represent the maximum values of h i and h ij , respectively, let us see how the intersection of these limits and the green region should be interpreted. We specifically focus on the baryogenesis line that is inside this region where the model produces the observed DM relic abundance and the desired baryon asymmetry. In both panels, the baryogenesis line is above the double proton decay bound and yet consistent with it. We only need h 1 h 12 to not be the largest combination of these couplings on this line, i.e., that h 1 and/or h 12 are sufficiently smaller than the largest values of h i and h ij respectively. As for the n −n oscillation limit, any τ n−n = const. contour in the orange band is compatible with all points in the green band that are above it. Again, it will be enough if h 1 and/or h 13 are sufficiently smaller than the largest values of the h i and h ij couplings respectively. On the portion of the baryogenesis line that is within the orange band, both h 1 and h 13 can take their largest values while remaining compatible with the current bound. This results in 4 × 10 9 s ≤ τ n−n ≤ 5 × 10 10 s (in the left panel) and 7 × 10 8 s ≤ τ n−n ≤ 5 × 10 10 s Finally, we comment on the prospect to experimentally probe DM in our model via direct and indirect searches. Couplings of χ to the quarks, see Eq. (1), result in both spin-independent and spin-dependent interactions between DM and nucleons. The spin-independent elastic scattering cross section is ∝ h 1 4 /m 8 X [17]. The resulting σ SI within the allowed parameter space in Fig. 1 is many orders of magnitude below the current bounds from direct detection experiments [35]. The spin-dependent cross section is found to be ∝ h 1 4 /m 4 X [17]. While σ SD σ SI , it is still much smaller than the current limits from direct detection searches [36]. Also, the DM annihilation rate at the present time follows σ ann v ∝ h 1 4 /m 4 X [17]. This is far too small to give rise to any detectable signals in indirect detection searches [37,38] even after enhancement of the small-scale DM perturbations during EMD is taken into account [8,10,12] 12 .

V. CONCLUSION
In this paper, we presented a minimal extension of the SM that gives rise to the observed DM relic abundance and accommodates successful baryogenesis. This model involves an O(TeV) colored scalar X and two Majorana fermions N and χ that are SM singlets. N is a neutral LLP with m N ∼ O(100) GeV and χ is the DM candidate whose stability requires that m χ ≈ m p . The correct DM abundance and baryon asymmetry are obtained through an interplay between thermal and nonthermal effects.
Starting in a RD Universe at T m X , decays of X bring N to a thermal abundance while also producing DM in the freeze-in regime. N subsequently dominates the energy density of the Universe as a result of its long lifetime and drives an epoch of EMD. N eventually decays and establishes a RD Universe before the onset of BBN. N decay is an additional source of DM production and, due to its baryon-number-violating nature, also gives rise to baryogenesis. We found the allowed parameter space where the observed DM relic abundance and baryon asymmetry are obtained while satisfying other phenomenological and cosmological constraints. This region of the parameter space overlaps with the MATHUSLA sweet spot for hadronically decaying LLPs, which warrants a careful study of the discovery prospects of our model. It can also be probed by next-generation n −n oscillation experiments, and a future positive signal can provide further information about the flavor structure of the baryon-number-violating couplings of the model.
The model we presented here is robust and largely independent of the postinflationary thermal history as it only requires that the Universe be in a RD phase at T O(TeV). This work can be extended in a number of directions. A natural question on the theoretical side is its embedding in UV-complete models of the early Universe like those arising from string theory. Our model can work in this set up, which typically predicts epochs of EMD driven by string moduli, as long as the last phase of modulus-driven EMD reheats the Universe to T m X . Models with high-scale supersymmetry breaking (for example, see [39]) can be a suitable framework in this regard. On the phenomenological side, it would be interesting to study dijet and monojet signals associated with X decay at the LHC (similar to the analysis performed in [23]) as a complementary probe of the model parameter space. We leave a detailed investigation along these lines for future work.
Here, f X and f q are the occupation numbers of X and q respectively. The factor of 3 on the RH side of these equations accounts for the fact that X is a color triplet. As shown in [13], these equations can be solved to find: where γ N ≡ Γ X→N /H(T = m X ) and γ χ ≡ Γ X→χ /H(T = m X ). Here, T i is the initial temperature at which we set f N ( p) = 0. We can then calculate the comoving number densities of N and χ, n co N (t) and n co χ (t) respectively, as functions of time: n co where the factor of 2 counts the internal degrees of freedom of N and χ and In this expression, γ without a subscript is understood to be either γ N or γ χ as it is the same for both.
In the left panel of Fig. 3, we show F (γ, T ) as a function of T /m X . Note that it saturates before T ≈ 0.1m X regardless of the value of γ. We therefore take F (γ) = F (γ, T 0.1m X ) to be the final comoving number density produced by X decay normalized to the equilibrium value. The right panel of Fig. 3 shows F (γ) as a function of γ, and we see that this function is roughly equal to F (γ) ≈ 7γ in the freeze-in regime (matching the determination in Ref. [40]), before flattening to F (γ) ≈ 1 once γ 1. Furthermore, there is no dependence on the initial time as long as it is before T ∼ m X .

FIG. 3:
Evolution of the comoving number density of either N or χ produced by X decay, normalized to the equilibrium number density. The left panel shows evolution with temperature for different values of γ while the right panel shows the dependence of F on γ at low temperatures (as compared to m X ) after it has saturated.

B. Contribution of Thermal Bath to DM Relic Abundance
Here, we show that the dominant contribution to the DM relic abundance from the thermal bath comes from X decay. As mentioned, we have F (γ χ ) 1, and hence DM production from annihilations in the thermal bath occurs in the freeze-in regime. X quanta are kinematically accessible to the thermal bath at T m X , and decay of on-shell X quanta (which we discussed in the text) dominates over X-mediated annihilations. At T dom T m X as well as in the adiabatic part of EMD 13 the temperature is simply redshifted according to T ∝ a −1 . The abundance of DM particles produced from annihilations in this interval follows (see [41]): DM production in the entropy-generating phase of EMD, where T ∝ a −3/8 , mainly occurs at T ∼ m χ /4 resulting in (for example, see [8]): The rates for DM production from tree-level qq/qq/qq annihilations (via X exchange) and qq/ll annihilations at 13 Radiation from the prior RD phase dominates over that produced from N decay in this period.
the one-loop level (via photon exchange) respectively follow: Using these expressions, we see that the RH sides of Eqs. (32,33) are totally overwhelmed by the first term on the RH side of Eq. (19) throughout the parameter space shown in Fig. 1. This confirms that X decay makes the main contribution of the thermal bath to the DM abundance in the allowed parameter space. One may also consider N q ↔ χq scatterings mediated by X that could in principle establish chemical equilibrium between N and χ. However, for the values of h in Fig. 1, the rate for this process is well below the Hubble rate at temperatures T < m X . This is expected since if X decay cannot bring χ into equilibrium, the scattering process, which has a smaller rate, will be even less efficient.

C. Constraints from DM Free Streaming
To see the importance of free streaming of DM particles, we first find the rate for scattering of DM particles off the SM particles in the thermal bath at the time of their production from N and X decay.
• DM from N decay. DM particles produced from N decay have a typical momentum p m N /2 when most production happens at H ∼ Γ N . The rates for χ scattering off u quarks (via X exchange at tree-level) and scattering off quarks and leptons (via photon exchange, which arises at the one-loop level) in the thermal bath respectively follow: where (m N T dec ) 1/2 in the numerator of the first expression is the momentum in the center-of-mass frame at H ∼ Γ N .
• X decay. Another source of DM production is X decay. As shown above, for F (γ χ ) 1, the comoving number density of χ is saturated at T 0.1m X . We can therefore take the typical momentum of DM particles at this time to be p ∼ m X /2. The rates for χ scattering off u quarks (via X exchange) and scattering off quarks and leptons (via photon exchange) at that time are respectively given by: The rate to maintain kinetic equilibrium Γ kin is much smaller than Γ scatt as many scatterings are needed for DM to exchange enough energy with the thermal bath (for example, see [42]). This implies that while Γ scatt > H is a necessary condition, it will not be sufficient to ensure that DM particles are in kinetic equilibrium with the thermal bath. On the other hand, DM particles are kinetically decoupled from the thermal bath whenever Γ dec < H.
From Eq. (36) we see that Γ scatt (T = T dec ) H dec , where H dec ∼ T 2 dec /M P , throughout the allowed parameter space shown in Fig. 1. This implies that DM particles produced from N decay are kinetically decoupled at H dec , and hence their typical momentum at H dec is p ∼ m N /2. Similarly, it turns out from Eq. (38) that Γ scatt (T = 0.1m X ) H(T = 0.1m X ) in the allowed parameter space. Thus, DM particles produced from X decay have an initial momentum p ∼ m X /2 that is redshifted ∝ a −1 due to the expansion of the Universe. As a result, their momentum is p T dec at the end of the EMD epoch when the Hubble rate is H dec 14 . A detailed analysis in [11] has constrained the boost factor of DM particles at the end of EMD as a function of T dec based on their impact on the matter power spectrum. For fixed values of m N and m X , this translates into a bound in the h − h plane. For the benchmark points discussed in Section III, X decay is the dominant source of DM production within the allowed parameter space, and hence the typical boost factor of DM particles at H dec is small. As a result, the constraint from the suppression of the matter power spectrum is easily satisfied throughout the entire parameter space shown in Fig. 1.
For heavier masses, however, limits from free-streaming may become relevant. In Fig. 4 we show the h − h plane for a high-mass case with m X = 300 TeV and m N = 10 TeV, as well as a dash-dot magenta line that corresponds to the 2σ confidence bound from [11] on the scale of suppression of the matter power spectrum, which is related to the free-streaming length. We see that in this high-mass case, our allowed parameter space is partially constrained from the bottom-left by these considerations. This is particularly relevant in the top-left corner where DM is produced primarily from N decay as can be seen from the turnover in the h contours. In the rest of the green region, most of the DM particles are produced from X decay and have a much smaller boost factor (∼ T dec /m χ instead of ∼ m N /m χ ).