Early Matter Domination from Long-Lived Particles in the Visible Sector

We show that a nonstandard cosmological history with a period of early matter domination driven by a sub-TeV visible-sector particle can arise rather naturally. This scenario involves a long-lived standard model singlet that acquires a thermal abundance at high temperatures from decays and inverse decays of a parent particle with SM charge(s), and subsequently dominates the energy density of the Universe as a frozen species. Entropy generation at the end of early matter domination dilutes the abundance of dangerous relics (such as gravitinos) by a factor as large as $10^4$. The scenario can accommodate the correct dark matter relic abundance for cases with $\langle \sigma_{\rm ann} v \rangle_{\rm f} \lessgtr 3 \times 10^{-26}$cm$^3$s$^{-1}$. More importantly, the allowed parameter space can be directly probed by proposed searches for neutral long-lived particles at the energy frontier, allowing us to use particle physics experiments to reconstruct the cosmological history just prior to big bang nucleosynthesis.


II. EARLY MATTER DOMINATION FROM THE VISIBLE SECTOR
Our scenario involves a SM singlet in the visible sector, N , that reaches equilibrium at early times via interactions with the thermal bath. It behaves as a matter component after becoming nonrelativistic and, provided that it is sufficiently long lived, dominates the energy density of the Universe. This results in an epoch of EMD that ends when N decay establishes a RD Universe and sets the stage for BBN.

A. The Scenario
For a quantitative study of this scenario, let us consider an extension of the SM with the following Lagrangian: Here ψ collectively denotes the SM fermions, X is a scalar with mass m X that is charged under the SM, and N is a Majorana fermion with mass m N m X that is a SM singlet. Since the SM is a chiral theory, gauge invariance implies that each ψ in Eq. (1) has a certain chirality. At energies E m X , which is relevant for N decay, the above Lagrangian gives rise to an effective four-fermion interaction hh † N ψψψ/m 2 X 1 . This results in N decay into three-body final states consisting of the SM fermions.
Starting in a RD phase at T m X , established at the end of inflationary reheating (for reviews, see [17,18]) or from the decay of a heavy modulus, the important stages of the thermal history in our scenario, arranged by decreasing Hubble rate, are as follows: (1) H H(T = m X ). At this stage, X is in thermal equilibrium due to its gauge interactions with the SM particles. N also acquires a thermal abundance, n N ∝ T 3 and ρ N ∝ T 4 , provided that the partial decay width of X to N satisfies: This ensures that X decay and inverse decay reach equilibrium before X becomes nonrelativistic (see Appendix A for details).
(2) H(T = m X ) > H H(T = m N ). During this stage, N particles are relativistic and n N ∝ a −3 (where a is the scale factor of the Universe). More precisely, g * n N ∝ a −3 as long as N is in kinetic equilibrium with the thermal bath, while n N ∝ a −3 when it is kinetically decoupled. Note, however, that the difference is negligible since the number of relativistic degrees of freedom, g * , changes minimally (if at all) for m N T < m X .
(3) H(T = m N ) > H H dom . N quanta become nonrelativistic during this stage. Their comoving energy density becomes frozen and remains constant provided that the rate for N N → ψψ and N ψ → ψψ satisfy: These relations ensure that N self-annihilation and its annihilation with SM particles are inefficient at T m N . Once the energy density of N becomes comparable to that of radiation, ρ N ρ R , N starts to dominate the energy density of the Universe. This happens at H = H dom , which can be found from the following: We have used the fact that ρ N ∝ a −3 , while ρ R ∝ g −1/3 * a −4 , and a(t) ∝ H −1/2 during RD 2 . The number of relativistic degrees of freedom at T = m N and T = T dom are given by g * N and g * dom respectively. The factor of 2 on the LH side 1 Examples of gauge-invariant effective interactions, in terms of left-handed (LH) chiral fermions, include N u c u c d c , N LLe c , and N QLd c .
These arise when X is a color triplet, an electroweak doublet, and a scalar leptoquark respectively. We will discuss the first case in detail shortly. 2 Recall that the comoving entropy density of the Universe is constant during RD. This implies that g * T 3 ∝ a −3 , and hence ρ R ∝ g * T 4 ∝ of Eq. (5) accounts for the two degrees of freedom associated with the Majorana fermion N . Therefore, we arrive at: The necessary condition for N dominance is: where Γ N is the width of N decay into three SM particles.
The Universe is in an EMD era driven by N during this stage, which ends when N quanta decay. N decay reheats the Universe to a temperature T dec . Note that N decay must complete before the onset of BBN, and hence T dec > T BBN 3 MeV 3 , which implies that: The decay releases entropy and dilutes any preexisting relic abundance by the following factor: Here, s before is the entropy density of radiation existing from prior stages at H = Γ N , while s after is the entropy density generated by N decay. They are given by: where we have used the fact that a(t) ∝ H −2/3 during the EMD epoch. The Universe enters a RD phase at H = Γ N , and hence T dec (90/π 2 g * dec ) 1/4 (Γ N M P ) 1/2 . We therefore find: B

. An Explicit Realization
We now present a specific realization of the scenario introduced above. It is based on a minimal extension of the SM that was proposed for low-scale baryogenesis and DM [20] (see [21] for a supersymmetric version). The Lagrangian (using two-component Weyl fermions) is: Here, u c and d c stand for the LH up-type and down-type antiquarks respectively. i and j are flavor indices, color indices are omitted for simplicity, and h ij is an antisymmetric tensor. X is an iso-singlet color-triplet scalar of hypercharge +4/3. N is a singlet fermion, with m N m X , which may be charged under a higher-ranked gauge group that includes the SM.
At energies E m X , one finds an effective four-fermion interaction N u c i d c j d c k after integrating out X. This results in N decay to three quark, and since N is a Majorana fermion, three antiquark final states. N decay can be the origin of baryon asymmetry in the model as described in [21]. Moreover, this model can explain the DM content of the Universe. In the supersymmetric version, the scalar partner of N is a natural DM candidate [21,22]. In the nonsupersymmetric version, a second copy of N that has approximately the same mass as the proton is a viable DM candidate [20]. This can in addition address the baryon-DM coincidence puzzle.
This model also has testable predictions for phenomenology. The baryon-number-violating couplings of X result in n −n oscillation. It also has novel monojet/monotop signals, as well as dijet events, at the LHC [23][24][25]. The most stringent experimental bound on the model parameters comes from double proton decay pp → K + K + . For m X ∼ 10 TeV and m N ∼ 1 TeV, so that the model will be testable at colliders, the resulting limit is |h 1 h 12 | < 10 −4 [26]. There are also bounds from K 0 s -K 0 s and B 0 s -B 0 s mixing, but they can be easily satisfied in the allowed parameter space.

III. RESULTS
In this section we present our results. First, we would like to demonstrate that all of the conditions needed for the success of our scenario, given in Eqs. (2,3,4,7,8), can indeed be simultaneously satisfied. Note that, as explained before, ψ particles are SM fermions of a specific chirality. For simplicity, we consider X couplings to given flavor combinations in N ψ and ψψ. The first term in the Lagrangian in Eq. (1) results in the following (rest frame) partial decay width for X (provided that m N m X ): The same term also leads to the following rate for N self-annihilation at energies E m X : Here, C 1 is a multiplicity factor that is associated with gauge charges of the final state ψψ * . Combination of the two terms in Eq. (1) allows N annihilation with ψ at the following rate (at energies below m X ): where C 2 is a multiplicity factor that counts all possible combinations of gauge charges in the initial and final states. The factor of 3 on the RH side is due to the fact that N can annihilate with any of the three ψ's that show up in the effective interaction N ψψψ. For the number densities in Eqs. (14,15) we use the equilibrium value for a relativistic fermionic degree of freedom 3ζ(3)T 3 /4π 2 . Finally, the width for three-body decay of N in its rest frame is: where the factor of 2 on the RH side arises because N , which is a Majorana fermion, decays to both ψψψ and ψ * ψ * ψ * final states.
To put things in perspective, consider the model described by Eq. (12). For simplicity, we consider X coupling to a single flavor combination u c i d c j d c k where i, j, k are fixed (with j = k). In this model, X → N u i decay brings N into equilibrium at T m X . The multiplicity factor for N self-annihilation in Eq. (14) is C 1 = 3 because of the color of u i in the final state. The multiplicity factor for N annihilation Eq. (15) and N decay Eq. (16) is C 2 = 6 due to possible color combinations in u i d j d k . In general, C 1 and C 2 are model-dependent factors that can be found for a given model. Also, we need to sum over all flavor combinations of ψ's that contribute to the various rates in Eqs. (13)(14)(15)(16).
We now impose the conditions given in Eqs. (2,3,4,7,8) by using the expressions in Eqs. (13)(14)(15)(16). As a proof of concept, we set C 1 = C 2 = 1 to demonstrate that all of these conditions can be simultaneously satisfied. In Fig. 1, we show the allowed region in the h − h plane for various values of m X and m N that are in the ballpark for testability of the scenario at colliders 4 . We see that the conditions in Eqs. (4,7), corresponding to the red and blue lines marked 2 and 1 respectively, essentially lie on top of each other. This can be understood upon substitution of Eqs. (6,15,16) as these conditions have exactly the same dependence on h, h , m X , and m N up to an overall numerical factor.
In Fig. 2, we show the allowed region of the m N − τ N plane for the values of m X used in Fig. 1, with the same coloring and labeling of the various conditions. We note that the condition in Eq. (2) has no dependence on m N or τ N , and when Eq. (2) is satisfied by even a few orders of magnitude, the condition in Eq. (3) is comfortably satisfied as well, as seen in Fig. 1. Therefore we do not show either of these conditions in Fig. 2 (instead, we discuss them in terms of the m N − h plane in Appendix C). The condition in Eq. (8) is readily translated into an upper bound on τ N ≡ Γ −1 N , and is shown by the yellow vertical line marked 3. The conditions marked 1 and 2, shown in blue and red, again lie nearly on top of each other for most of the parameter space. However, they now deviate from straight lines in the bottom-right of both panels of Fig. 2 due to changes in g * (T ) (see Appendix B). In particular, as seen in Eq. (6), the blue curve depends on both g * N and g * dom which, along with the sharp drop in g * (T ) near the QCD phase transition, is responsible for the behavior at low m N . One comment is in order at this point. When m N is very close to m X , the expression in Eq. (13) is no longer a good approximation, and Γ X→N decreases due to phase-space suppression. Moreover, as shown in Appendix A, the comoving energy density of N reaches a thermal value and then freezes at a temperature T m X /5 provided that m N m X , where the exact temperature depends on the ratio Γ X→N /H(T = m X ). Together, these considerations imply that as m N approaches m X , the final abundance of N and hence the onset of the EMD epoch will change as compared to the m N m X limit. For this reason, the portion of the allowed region in Fig. 2 (as well as in subsequent figures) that is very close to the m N = m X line may be considered "less viable".

IV. CONNECTION TO EXPERIMENT
In this section, we discuss how the scenario may be tested experimentally. Specifically, we explore possible correlations between the allowed regions of the m N − τ N plane and σ ann v f in light of DM indirect detection searches. We also point out the prospects for detecting N via proposed searches for long-lived particles at the LHC.  (4,7,8), as well as m N < m X . The various conditions are labeled as in Fig. 1, with the addition of condition 6. The faded portion at the top of the allowed region, just below the m N = m X line, emphasizes that m N m X is needed for our scenario.

A. Implications for Dark Matter
If DM freeze-out occurs after the end of the EMD epoch driven by N , then T f < T dec and the standard picture of thermal DM will hold. Otherwise, EMD will affect the DM abundance. If freeze-out happens before the onset of EMD, T f > T dom , the relic abundance will be set during RD prior to EMD, and entropy release by N decay will subsequently dilute it. If freeze-out occurs during EMD, T dec < T f < T dom , thermal production of DM will also be altered [9,27]. Moreover, DM particles can be directly produced in the decay of the matter-like component that drives EMD [28][29][30][31].
In our scenario, the EMD epoch starts quite late and can end just before the onset of BBN. Eq. (6) results in T dom ∼ m N /50, where we have used g * ∼ O(100), whereas T dec 3 MeV. On the other hand, T f ∼ m DM /20. Therefore, unless DM is very light, its relic abundance will be affected by EMD. Below, we discuss the cases with small and large σ ann v f separately.
In a standard thermal history, this leads to overproduction of DM. However, an epoch of EMD can regulate the relic abundance. Let us consider the case where m N < m DM . As mentioned above, this implies that T dom < T f , and hence freeze-out occurs in the RD phase preceding the EMD era. Since N is lighter than the DM particle, direct production of DM from N decay is kinematically forbidden. Therefore, the only role of EMD is to dilute the overabundance of DM. The final DM abundance, normalized by the entropy density s, then follows: where the dilution factor d is given by the expression in Eq. (11), and the observed DM relic abundance is: Thus, to obtain the correct relic abundance, we need to have: From Eqs. (6,11), we see that d is basically determined once we know m N and Γ N . Using this, in Fig. 3 we plot contours of σ ann v f in the m N − τ N plane that yield the observed DM abundance for the same values of m X as in Fig. 2. The contours show σ ann v f decreasing by factors of 10, relative to the nominal value of 3 × 10 −26 cm 3 s −1 , corresponding to successively longer periods of EMD with more and more dilution. The region of interest is where the dashed contours overlap with the green allowed region of our scenario. We see that the scenario can give the desired relic abundance for σ ann v f as small as ∼ 10 −30 cm 3 s −1 , which is well below the current Fermi-LAT bounds from observations of dwarf spheroidal galaxies [32] and newly discovered Milky Way satellites [33].  Fig. 2. This essentially corresponds to the absence of EMD as Γ N = H dom in this case, thus recovering the standard freeze-out scenario. As before, the temperature dependence of g * (T ) is responsible for deviations from straight lines.
2. σannv f > 3 × 10 −26 cm 3 s −1 Freeze-out in a standard thermal history, as well as in an EMD scenario, leads to an underabundance of DM in this case. However, direct production of DM from N decay combined with a large value of σ ann v f can yield the correct relic abundance. If N decay results in a DM density that is higher than the observed value, and σ ann v f is sufficiently large, residual annihilation at the end of the EMD epoch (when H Γ N ) gives [28,29,[34][35][36]: This matches the observed DM abundance provided that: We note that σ ann v f > 3 × 10 −26 cm 3 s −1 is allowed within the 20 GeV m DM 100 TeV mass range. For DM masses up to a few TeV, the upper limit on σ ann v f is set by a recent analysis [37] of Fermi results (unless there is P -wave annihilation or coannihilation), while for larger values of DM mass, it is set by the well-known unitarity bound [38]. The largest allowed value is σ ann v f 2 × 10 −23 cm 3 s −1 , which happens at the intersection of these constraints. Fig. 4 depicts regions in the m N − τ N plane that yield the correct relic abundance for a given m DM . We use the allowed range for σ ann v f from combined experimental and unitarity bounds in [37] for a given m DM , and translate it into a range for T dec using Eq. (21) with T f ∼ m DM /20. T dec is readily related to Γ N which then restricts τ N to be within a vertical band. By imposing m N > m DM , so that N decay to DM is kinematically allowed, we find the allowed regions shown in purple. Different panels in Fig. 4  region of interest is where the vertical band overlaps with the allowed region of our scenario. The only dependence on m X is in the orange horizontal line corresponding to m N = m X . As m X increases, this line moves upward resulting in more overlap between the two shaded regions.

Other Considerations
The EMD epoch driven by N can also regulate the DM relic abundance in cases that DM production is not associated with processes in the thermal bath. An important example is axion DM in the form of coherent oscillations of the QCD axion arising due to an initial displacement from its minimum. EMD dilutes axionic DM if T dec < Λ QCD 200 MeV thereby avoiding any overproduction [39]. The maximum dilution is obtained for T dec T BBN 3 MeV, which allows an axion decay constant of order f a 10 14 GeV without fine tuning of the initial misalignment angle. In supersymmetric extensions of the SM, this can be achieved with dilution from the decay of the axino and/or saxion [40,41], however, our scenario works just as well in a nonsupersymmetric set up.
Our scenario can be particularly helpful with DM at the high end of the mass spectrum, i.e., m DM O(TeV). In this case, according to Eq. (18), obtaining the observed relic abundance requires a very small value of n DM /s, which could pose a challenge to any production mechanism of superheavy DM. However, entropy release by N decay can dilute the existing DM density by a few orders of magnitude thereby enlarging the parameter space of any model with extremely heavy DM candidates (for example, see [42]).
Moreover, the EMD epoch driven by N dilutes any relic abundance produced at prior stages of the cosmological history. A notable example is the observed matter-antimatter asymmetry of the Universe. Some mechanisms, notably Affleck-Dine baryogenesis [43,44], can generate an excessively large baryon asymmetry. Entropy released in N decay can help regulate the baryon asymmetry in such cases. The same applies to dangerous relics such as unstable gravitinos of weak-scale mass whose abundance is constrained by the success of BBN [45,46].

B. Testable Collider Signals
Due to its long lifetime, required for the success of our scenario, N is an example of a long-lived particle (LLP). This provides us with an opportunity to directly probe a sub-TeV N via the proposed LLP searches at the LHC. As we have seen, the allowed region in the m N − τ N plane covers a wide range of N lifetime. The rest-frame lifetime τ N of N particles produced from X decay in colliders is related to the decay length l N according to l N =bcτ N . Here,b is the average boost factor of N , which we take to beb ∼ m X /2m N . Decay lengths above 10 −2 cm are within the reach of the LHC. For 10 −2 cm < l N < 10 2 cm, N production and decay gives rise to displaced vertices at the LHC, while the case with 10 2 cm < l N < 10 4 cm leads to displaced jet/lepton signals.
However, neutral LLPs (like our N ) with decay lengths above 100 m are particularly difficult to probe because of the limited sensitivity of the LHC main detectors. The recently proposed MATHUSLA (MAssive Timing Hodoscope for Ultra Stable neutraL pArticles) detector concept [14] is a minimally instrumented, large-volume surface detector located near ATLAS or CMS. It would search for neutral LLPs produced in the high luminosity LHC (HL-LHC) collisions, extending the lifetime range by a few orders of magnitude compared to the main detectors. It could discover LLPs with decay lengths up to 3 × 10 7 m, which correspond to lifetimes close to the age of the Universe at the onset of BBN [15].
The case when l N > 100 m is particularly interesting as it may be probed by ATLAS and CMS as well as MATHUSLA. To be specific, let us focus on the model in Eq. (12) with m X ∼ 1 − 10 TeV and m N ∼ 100 GeV − 1 TeV as a canonical example. X can be produced at the LHC due to its couplings to the quarks in this model, and its decay yields dijet and monojet signals (the latter accompanied by N ). These signals have been studied in detail in [23,25], and in a recent analysis by CMS [47], when N is absolutely stable and m N ≈ 1 GeV (hence light DM). One may generalize these studies to m N ∼ 100 GeV − 1 TeV in a rather straightforward manner.
Interestingly, our benchmark point also overlaps with MATHUSLA's most important physics target, namely hadronically decaying LLPs with masses in the O(10 GeV) − O(100 GeV) range [16]. The collaboration has studied the discovery potential of such LLPs produced from exotic Higgs decays [15], or through mixing of the (scalar) LLP with Higgs [16]. In our model, as mentioned, the main production channel for N is decay of X particles. This should be implemented properly to calculate the cross-section for N production. Thus, combining MATHUSLA with the main detectors will allow us to probe the m N − τ N plane at the HL-LHC. Fig. 5 shows the allowed region of our scenario in the m N − τ N plane with regimes of the decay length and corresponding prospects for detection at the LHC. Though part of our allowed region extends to l N < 10 2 cm, most of it lies at longer lifetimes with l N > 10 4 cm. We also see that the region relevant for MATHUSLA has some overlap with the regions that yield the correct DM abundance for σ ann v f > 3 × 10 −26 cm 3 s −1 in Fig. 4. The situation is even better, see Fig. 3, in the case that σ ann v f < 3 × 10 −26 cm 3 s −1 . Note that we have expressed the average boost factor of N particles asb ∼ m X /2m N , which is strictly correct when X decays at rest. However, X particles produced at the LHC are boosted themselves, leading to a largerb for N . In this case, the lines corresponding to different values of l N would move up and to the left, making the region within the reach of MATHUSLA even larger.
Potential collider measurement of m N and τ N has very interesting cosmological implications. τ N (equivalently Γ N ) is a direct measure of the temperature T dec when N decay reheats the Universe. This leads to the tantalizing possibility of determining the highest temperature of the Universe in the last phase of RD, which is relevant for BBN, via particle physics experiments. By knowing m N , we can also find the onset of the EMD period that is driven by N through Eq. (6). Therefore, in principle, we can directly probe the entire EMD epoch ending just before the onset of BBN with collider experiments. The diagonal band between 10 −2 cm < l N < 10 2 cm leads to displaced vertices at the LHC. The next band, with 10 2 cm < l N < 10 4 cm, leads to displaced jet/lepton signals. Finally, the right-most region, with l N > 10 4 cm, would be probed by MATHUSLA. The darker shaded band within the allowed region of our scenario corresponds to m N = 10 GeV−100 GeV, which is the most important MATHUSLA target for hadronically decaying LLPs.

V. CONCLUSION
In this paper, we presented a nonstandard cosmological history scenario where a visible-sector particle drives a period of EMD. This scenario involves a sub-TeV SM singlet N that acquires a thermal abundance at temperatures well above its mass as a result of the decay/inverse decay of a parent particle X with SM charge(s). Then, being long lived, N dominates the energy density of the Universe as a frozen species and leads to an EMD epoch that can last until the onset of BBN. The scenario works for fermionic and bosonic N equally well. Moreover, since the energy density of N originates from the thermal bath, it automatically evades isocurvature bounds. We discussed an explicit realization of the scenario where N is a Majorana fermion coupled to quarks.
We outlined the necessary conditions in order for N to reach and maintain an equilibrium energy density, dominate the Universe, and decay before the onset of BBN. We showed that these conditions can be simultaneously satisfied in significant parts of the parameter space. It is possible to obtain the correct DM relic abundance within this allowed parameter space for both cases with σ ann v f < 3 × 10 −26 cm 3 s −1 and σ ann v f > 3 × 10 −26 cm 3 s −1 . Moreover, the entropy release from N decay can regulate baryon asymmetry generated at earlier stages of the cosmological history and dilute the abundance of dangerous relics like unstable gravitinos.
An interesting aspect of this scenario is that the m N − τ N plane may be directly probed by collider experiments. In large parts of the allowed parameter space, τ N is in the range for discovery by the MATHUSLA proposal at the HL-LHC. This is particularly the case for a hadronically decaying N with m N ∼ O(100 GeV), which is MATHUSLA's most important physics target. Measuring τ N and m N exhibits a profound interplay between particle physics and cosmology. Determining τ N will give us the highest temperature in the RD phase that sets the stage for BBN. Using information about m N , in tandem, we can fully reconstruct the EMD era driven by N .
In summary, a visible sector particle N with sub-TeV mass can give rise to a (final) period of EMD in the early Universe rather naturally. This scenario is robust, and largely independent from the details of the earlier stages of the cosmological history, as long as the Universe is in a RD phase at T m N . A logical next step is embedding the scenario in realistic extensions of the SM with TeV-scale long-lived fermions or scalars. Another important direction is to employ exciting rapid developments in LLP searches and perform detailed analysis of the discovery prospect of such models at the energy frontier. We leave these investigations for future work.

A. Evolution of the Energy Density
Here, we would like to derive the time evolution of the energy density of N particles produced via the interaction term hXN ψ in Eq. (1). We work in the limit m ψ m N m X . For simplicity, we take m N = m ψ = 0 below. This is a good approximation as long as we are interested in the evolution of ρ N over time scales where T m N . Let us start with the equation that governs the occupation number of N , denoted by f N : Here, f X and f ψ are the occupation numbers of X and ψ respectively, and M is the Feynman amplitude for the X → N ψ * process (note that, being a Majorana fermion, N is its own antiparticle). Since X and ψ are in thermal equilibrium, f X = f eq X and f ψ = f eq ψ . We also have |M| = hm X . Integration over p ψ yields: where: Conservation of momentum implies that: Using spherical coordinates, with the z axis chosen in the direction of p N , we have: where p X ≡ | p X | and p N ≡ | p N |. Integrating over the solid angle, and noting isotropy about p N , we find: where Γ X→N = h 2 m X /16π, see Eq. (13). We have used f eq ψ = f eq X /f eq N , with f eq X ( p X ) = e −βE X and f eq N ( p N ) = e −βp N . We note that p X dp X = E X dE X , and hence: Here, E X,min is the minimum value of E X that results in a given value p N due to X decay. Conservation of energy and momentum together, see Eq. (26), imply: This gives: Thus, the minimum of E X occurs when θ = 0 or θ = π. The former corresponds to forward production of N and is the case when p N > m X /2. The latter corresponds to backward production of N , which happens to be the case when p N < m X /2. In both cases, we find: We see that, as expected, E X,min ≥ m X . The lower bound is saturated when p N = m X /2, which occurs when the decaying X is at rest. With the final integration, Eq. (28) then takes the following form: Choosing an initial time t i when T = T i , we have p N (t) = p N (t i )a(t i )/a(t) and T (t) = T i a(t i )/a(t), where a is the scale factor of Universe with a(t) ∝ t 1/2 during RD. For simplicity, we neglect any changes in g * (T ) due to the high temperatures we are considering here. Accounting for expansion of the Universe, we thus have: This can be rewritten as: where f eq N ( p N ) = e −p N has no dependence on time. This allows us to perform the integral of both sides exactly. Starting with f N ( p N ) = 0 at t = t i and H i = 1/2t i , we then find: where γ ≡ Γ X→N /H(T = m X ) and we have made use of H(T = m X )/H i = α 2 .
We can now calculate the comoving energy density of N as a function of time, ρ co N (t), as follows: where the factor of 2 counts the internal degrees of freedom of N . Substituting the expression in Eq. (35) for f N ( p N ), and noting that p N ∝ a −1 , we arrive at the following relation: We note that the ratio of the comoving energy density of N to its equilibrium value is mainly sensitive to γ. Dependence on α only shows up through the lower limit of the integral overt that is related to the initial time.
In Fig. 6 we show the evolution of Eq. (37) with respect to the temperature of the Universe. The horizontal axis corresponds to T /m X because Eq. (37) depends on this ratio rather than T itself. In the left panel, we see that even for the smallest value allowed in Eq. (2), expressed as γ = 1, the comoving energy density of N reaches ∼ 0.8ρ eq,co N around the time when T ∼ m X /5. Thus in our scenario N acquires a thermal energy density at T m N as long as m N m X /5. However, this would not be the case for γ 1. The reason being that in this case the maximum comoving energy density established by X decay (and inverse decay) is ρ eq,co N as the X number density is Boltzmann suppressed. In the right panel, we see that the value of T /m X for which ρ co N ρ eq,co N is essentially insensitive to α, as pointed out above. The only effect of α is to determine the starting point of N production in our calculation.
In the above derivation, we have only considered N production from the decay of X (and the accompanying inverse decay). At temperatures T m X , inelastic scattering of N off SM particles N ψ ↔ ψ * ψ * might be in equilibrium too. This would be an additional source of N production from the thermal bath and thereby help ρ co N reach its equilibrium value even faster. Elastic scattering of N off SM particles N ψ → N ψ would only redistribute the energy among N particles, which could affect the evolution of individual f N ( p N ), but would not change the total energy density. In Fig. 2 we showed the four conditions labeled 1, 2, 3, and 6, however conditions 4 and 5 of Fig. 1 were not shown. Here we present a short discussion to better understand how conditions 4 and 5 affect the interplay between the h − h and m N − τ N planes, and to aid in interpreting Figs. 1 and 2 together.
We define γ ≡ Γ X→N /H(T = m X ), as before, and η ≡ Γ N N →ψψ * /H(T = m N ) such that η = 1 and γ = 1 correspond to conditions 4 and 5 of Fig. 1 respectively. The main issue is that a single value of m N on Fig. 2 corresponds to a range of h, γ, and η. Fixing m X = 10 TeV and m N = 1 TeV as an example, we can obtain the full range of h, γ, and η that are allowed for this single combination of the masses while satisfying both Eqs. (2,3). If we choose γ = 1 such that we sit on the line for condition 5 in the bottom-right panel of Fig. 1, we get h ≈ 10 −6 which then gives η ≈ 10 −16 . On the other hand, if we first choose η = 1 so that we are on the line for condition 4, we then have h ≈ 10 −2 and therefore γ ≈ 10 8 . The possible range of these three parameters for this single combination of m X and m N is therefore h ≈ (10 −6 − 10 −2 ), γ ≈ (1 − 10 8 ), and η ≈ (10 −16 − 1). In Fig. 8 we show conditions 4 and 5, as well as m N = m X , in the m N − h plane with the corresponding allowed region. Note that m N is restricted from below by the other conditions of our scenario, not shown here. For a given value of m N in Fig. 8, one can read off the allowed range of h and then calculate the corresponding values of the parameters γ and η using Eqs. (13,14).
We further note that conditions 4 and 5 meet at high m N , allowing γ = η = 1 to be simultaneously satisfied. However, this occurs far above m N = m X and is therefore out of reach unless m X is itself very large. Additionally, this intersection corresponds to a single value of h such that the entire allowed region of Fig. 1 becomes compressed to a line, as condition 4 moves down to join condition 5.