Early kinetic decoupling effect on the forbidden dark matter annihilations into standard model particles

The early kinetic decoupling (eKD) effect is an inevitable ingredient in calculating the relic density of dark matter (DM) for various well-motivated scenarios. It appears naturally in forbidden dark matter annihilation, the main focus of this work, which contains fermionic DM and a light singlet scalar that connects the DM and standard model (SM) leptons. The strong suppression of the scattering between DM and SM particles happens quite early in the DM depletion history, where the DM temperature drops away from the thermal equilibrium, $T_\chi<T_{\rm SM}$, leading to the decreased kinetic energy of DM. The forbidden annihilation thus becomes inefficient since small kinetic energy cannot help exceed the annihilation threshold, naturally leading to a larger abundance. To show the eKD discrepancy, we numerically solve the coupled Boltzmann equations that govern the evolution of DM number density and temperature. It is found that eKD significantly affects the DM abundance, resulting in almost an order of magnitude higher than that by the traditional calculation. We also discuss the constraints from experimental searches on the model parameters, where the viable parameter space shrinks when considering the eKD effect.


I. INTRODUCTION
Relic density is an essential topic for dark matter (DM) physics. The classical scenario that explains the observed abundance in the present Universe is thermal particle production in the early Universe, which is the so-called freezeout mechanism for the weakly interacting massive particle (WIMP) [1][2][3][4][5]. The DM particles were initially in thermal equilibrium with the heat bath via the intense interactions among them. The number density dilutes along with the expansion of the Universe and finally freezes out of the heat bath once the annihilation rate falls behind the cosmic expansion rate, resulting in a comoving constant.
Gondolo et al. [6,7] have developed the renowned treatment of calculating the DM relic density by solving the Boltzmann equation of the number density with high accuracy, called the "standard" method. One hypothesis entering this treatment is that DM keeps in local kinetic equilibrium with the thermal plasma during or even after the freeze-out process. The scatterings with the standard model (SM) particles have been at a much more intense level [8].
However, this is not always true for many well-motivated mechanisms, where kinetic equilibrium might decouple earlier than assumed, leading to the early kinetic decoupling (eKD) around the freeze-out period. The eKD effect has been extensively studied in the literature [9][10][11][12][13][14][15][16][17][18][19][20]. The most influenced scenarios include 1) resonant annihilation of dark matter [10,13,14], 2) Sommerfeld-enhanced annihilation [14], and 3) subthreshold annihilation (also known as the forbidden annihilation) [14]. In these regimes, the eKD occurs because the elastic scattering processes are suppressed; DM and the SM particles experience different temperatures. All these cases demonstrated the actual DM abundance can be affected * Electronic address: liuyu2020@s.ytu.edu.cn † Electronic address: xuewenliu@ytu.edu.cn ‡ Electronic address: zhubin@mail.nankai.edu.cn by up to 1 order of magnitude compared to the traditional method, at least in some parts of the parameter space.
Here, we study the eKD in the forbidden DM scenario [21]. In such a scenario, DM dominantly annihilates into heavier final states, which can proceed at finite temperatures in the early Universe, relying on the thermal tail with the high velocity of DM. Many works were devoted to studying the forbidden annihilations with different theoretical models and a variety of phenomenological topics [22][23][24][25][26][27][28]. In this work, we employ the model where forbidden DM annihilations into the SM leptons are mediated by a singlet scalar. Such channels were studied in Ref. [23], which is experimentally viable and predicts a very narrow mass range for DM, that can be tested at future beam dump experiments. However, the eKD effect has not been studied in the context of this model.
Actually, the scatterings of DM against SM particles are strongly suppressed in the early times because of the mass splitting between the DM and SM leptons. The DM temperature drops away from the thermal bath temperature, T χ < T SM , which leads to the decreased kinetic energy of DM. The forbidden annihilations thus become inefficient since small kinetic energy cannot help exceed the annihilation threshold. So this naturally leads to a larger abundance and the eKD effects in such a model should not be neglected. On the contrary, eKD will cause significant impacts on DM abundance. We investigate the relic density beyond the standard treatment used in Ref. [23] by considering the coupled Boltzmann differential equations, where the temperature evolution of the dark sector could be taken into account. We use the public code DRAKE [14] to perform the numerical calculations. We find a DM relic density that differs by up to an order of magnitude from the standard treatment and leaves a reduced feasible parameter space under the various experimental constraints.
The rest of this paper is structured as follows. In Sec. II, we start with a general description of the coupled Boltzmann equations that govern the evolution of DM number density and temperature. In addition, we discuss the DM model for forbidden annihilations and analyze the occurrence of early kinetic decoupling. Section III is devoted to a thorough study of DM relic density for the forbidden channels, and makes a detailed comparison between the new treatment and the traditional one. We further discuss various constraints from collider searches and astrophysical observations on the parameter space in Sec. IV. Finally, we conclude in Sec. V.

II. EARLY KINETIC DECOUPLING EFFECTS ON FORBIDDEN ANNIHILATIONS
A. Basic formulas Keeping kinetic equilibrium during and even after the freeze-out epoch is one underlying assumption for traditional relic density calculations. However, this is not always the case for various scenarios, where kinetic decoupling happens earlier than the chemical decoupling process. To study the DM relic density by taking into account the early kinetic decoupling effect, we should consider the following Boltzmann equation for DM phase-space distribution [10,14] where E is the energy of the DM, H is the Hubble constant, ⃗ p is the momentum of DM, and f χ is the DM phase-space density. The collision term C ann. represents the annihilation of DM particles into thermal bath particles, and C el. is for elastic scattering processes between DM and SM scattering partners. For two-body processes, where B and B ′ stand for particles in the thermal bath such as SM leptons, g χ is the number of internal degrees of freedom of DM, and f eq B is given by the Fermi-Dirac or Bose-Einstein distribution depending on the spin of B. The summation should be taken for all the internal degrees of freedom for all the particles. For the nonrelativistic DM, C el. can be simplified as the Fokker-Planck operator [8,[29][30][31] 1 : In the above, the momentum transfer rate γ(T ) is given by (see also Ref. [32]) where the differential cross section can be expressed as , and k 2 cm is given by Here E k is the energy of heat bath particle B. Note that k 2 During the chemical decoupling, the scattering processes may not be frequent enough to maintain the kinetic equilibrium, which means that DM particles own a different temperature T χ from the thermal plasma in their following evolution. A common definition of the DM temperature is which is also a function of the thermal bath temperature T . In this definition n χ is the number density of the DM, and s is the entropy density. Here y is a dimensionless version in analogy to the DM yield Y(= n χ /s).
To reach a suitable description of the DM temperature evolution and then explore the eKD effect on the chemical decoupling process, we should consider the second moment of f χ as a dynamical degree of freedom. By integrating Eq. (1) with g χ d 3 p (2π) 3 1 E and g χ d 3 p (2π) 3 1 E ⃗ p 2 E 2 , one obtains the zeroth and second moments of the Boltzmann equation, respectively. This leads to a relatively simple coupled system of Boltzmann differential equations ( denoted as the cBE method hereafter), Tχ . Note that the elastic scattering term given in Eq. (4) does not contribute to the zeroth moment term. This is a natural consequence because the elastic scattering processes do not change the number density of DM.
The above compact form of the differential equations contains the following thermally averaged cross sections, The thermal average ⟨σv⟩ 2,T is a variant of the commonly used thermal average ⟨σv⟩ T , and is explicitly stated in Ref. [10] and introduced as For ⟨σv⟩ T and ⟨σv⟩ 2,T , replace T χ by T in ⟨σv⟩ T χ and ⟨σv⟩ 2,T χ , respectively. And n eq χ (T χ ) is given by Tχ . (13) In this work, we use the numerical routine DRAKE to solve the coupled Boltzmann equations. The measured value of Ωh 2 by the Planck Collaboration is Ωh 2 = 0.120 ± 0.001 [33]. The viable parameter space is determined by matching this value.

B. Model and discussion on forbidden channels
We have adopted a simple model that only takes into account DM annihilations into SM leptons. The DM is a Dirac fermion coupled to the SM sector via the scalar portal ϕ. After the electroweak symmetry breaking, the effective Lagrangian can be written as where the indices on the couplings i, j = e, µ, τ. This model has been studied thoroughly for the forbidden mechanism (for detail, refer to Ref. [23]). The merits include the following: 1) DM mass is limited in quite a small window close to the masses of the SM leptons, which is a strong prediction that can be tested soon by colliders or beam-dump experiments. 2) Kinematically forbidden DM naturally evade the stringent constraints from the energy injections into the cosmic microwave background (CMB) [21]. In the forbidden DM scenario, the energy injection processes suffer the Boltzmann suppression at T ≲ eV, so that sub-GeV thermal relics are consistent with the experiment, making annihilations to SM leptons with DM masses m χ ≪ 10 GeV still viable. In this scenario, the DM relic density should be carefully scrutinized, as the eKD effect appears generic. The actual relic density receives a significant correction compared with the conventional method, as shown in the following sections.
Following Refs. [21,23], we consider a pair of DM particles dominantly annihilating into two SM particles 2ℓ with mass m χ < m ℓ . Cosmological constraints make forbidden annihilations into electrons unfeasible, including big bang nucleosynthesis (BBN) and CMB [34]. So, for simplicity, we only consider the µ + µ − and τ + τ − channels, with abbreviated couplings g (A) µ , g (A) τ , and g (A) χ , which allows us to explore all the relevant DM phenomenology systematically.
Applying the detailed balance condition for the DM number-changing process, the cross section of the forbidden channels is exponentially suppressed, where ∆ ≡ (m ℓ − m χ )/m χ . σ χ ≡ σ(χχ → ℓℓ ) , while σ ℓ is the cross section for the inverse process. When the annihilation rate becomes slower than the Hubble expansion, DM is no longer in equilibrium with the SM thermal bath, resulting in chemical decoupling.
What about the scattering between DM and SM particles during this period? From Eq. (5), the momentum transfer rate γ is proportional to an exponential factor: The full expression of the rates is listed in Appendix A. In the forbidden scenario, ∆ + 1 > 1 implies the scattering frequency experiences a strong suppression at a much earlier period. It is known that DM kinetically decouples out of the SM thermal bath as long as the momentum transfer rate γ is smaller than the Hubble expansion, γ < H. For illustration, we show the momentum transfer rate γ for the χµ ± → χµ ± scattering process, in Fig. 1. In comparison, we plot the Hubble parameter H(x) as a function of x. The red lines stand for the evolution of γ(x) and H(x) in the forbidden DM case where we take m χ = 0.1GeV. To demonstrate the distinctiveness of the forbidden DM, we also provide the results of a nonforbidden case where DM mass is larger than that of the annihilation products (m χ =1 GeV). The most remarkable finding is that γ of the forbidden case becoming comparable with H(x) happens much earlier than that of the nonforbidden case. The kinetic decoupling starts at around x = 20, which is usually the same time as DM chemically decoupled from the thermal bath. The reason for this very early kinetic decoupling is straightforward to understand as the result of an exponential suppressed momentum transfer rate, as derived in Eq. (16). One can obtain similar results for the χχ → τ + τ − forbidden channel. We conclude that eKD exists in the forbidden DM scenario. The next step is to find the effect of the eKD more concretely. We systematically study DM relic density in both cBE and traditional methods (denoted as nBE as in Ref. [14]) and then discuss the phenomenological possibilities.

III. RELIC DENSITY: COMPARISON BETWEEN CBE AND NBE APPROACHES
In this section, we compute the relic density in both the standard method (nBE) and the cBE approach. We restrict our study to the same range of m χ that 0.9m µ ≲ m χ ≲ m µ as derived in Ref. [23], in which forbidden annihilation into µ + µ − is experimentally viable. For the χχ → τ + τ − case, the corresponding DM mass is 0.8m τ ≲ m χ ≲ m τ . The detailed annihilation cross sections and the momentum transfer rates for scatterings are presented in Appendix A.
To reveal the effects of the eKD and the differences between the cBE and nBE approaches, we first find several benchmark points for the two forbidden channels, shown in Table I. We also set g χ = 0, (g A χ ) 2 /4π = 0.1 as in Ref. [23], for making a rough but straight comparison. The values of the rest model parameters are fixed to obtain the observed DM relic density for Ω cBE h 2 . With these inputs, we can compute the relic density in the nBE method. It can be seen that the ratio of Ω cBE /Ω nBE is sizable, even reaching an order of magnitude.
The significant difference between the cBE and nBE results exactly comes from the eKD effects. In Fig. 2, we show the temperature and abundance evolution for selected benchmarks in Tab. I. From the left panel, the green lines are the evolution curves of y χ , namely the temperature of DM, which depart from the thermal bath temperature (the gray line) at around x = 20. Qualitatively, DM needs higher momenta to overcome the annihilation threshold, leading to a self-cooling phase as soon as it is no longer kinetically coupled to the muons. This is why there is a drop and the temperatures evolve separately for the dark sector and the SM sector since then. DM annihilation becomes less efficient much earlier just because of this cooling, which results in a higher DM abundance than in the nBE approach, as shown in the right panel of Fig. 2. Note that the same cooling phenomena also have been found in Ref. [14]. In Fig. 3, we show a global picture of the eKD effect for the forbidden cases of χχ → µ + µ − and χχ → τ + τ − , where we define the deviations of the relic density in cBE and nBE approaches by Deviation ≡ (Ω cBE − Ω nBE )/Ω cBE . We display the results in the (m ϕ , g µ/τ ) plane by the density plotting method, with setting m χ = 0.1GeV and m χ = 1.48GeV for illustration. The sizeable deviations appear in most parameter spaces from 20% up to almost 100%. The maximum deviations seem to emerge in the resonance region. However, it should be noted that eKD effects already exist in the resonant annihilations of DM [10,13,14]. So in this region, one should study the eKD for resonance and forbidden channels together. It is beyond the scope of our work, as we mainly focus on the forbidden annihilations.
To emphasize the importance of the improved treatment of the decoupling history near the mass threshold, we plot in Fig.  4 the ratio of the resulting relic density to that of the standard nBE approach. Here the parameters satisfy the requirement for Ω nBE h 2 = 0.12. The different choices of coupling g A l correspond to different curves, as labeled in the plots. The mass ratio m ℓ /m χ can be divided into three regions. At lower mass ratios, dark matter evolves as ordinary WIMPs where kinetic equilibrium is maintained. The relic density derived by the cBE and nBE methods agrees with each other. The gray-shaded region is known as the resonance region, where 2m χ ≃ m ϕ . Here, the eKD effect arises due to the distinct cooling and heating effects of dark matter.For a more detailed discussion of the origin of these features, please refer to Ref. [10].
In the forbidden DM mass region (brown-shaded), we can see that the cBE results are larger than nBE several times with the same parameters. For µ + µ − case, 3 ≲ Ω cBE Ω nBE ≲ 10, and 2 ≲ Ω cBE Ω nBE ≲ 15 for the τ + τ − forbidden case. As the mass ratio m µ/τ /m χ increases, the coupling (g µ/τ ), which is needed to obtain the correct relic density, rises  FIG. 4: Relic density comparison between Ω cBE h 2 and Ω nBE h 2 for the χχ → µ + µ − and χχ → τ + τ − channels. Here, we choose the parameters to match the observed abundance in the nBE approach. The left diagram corresponds to the µ + µ − mode, and we take m ϕ = 0.26 GeV, and the right diagram presents the τ + τ − mode, for which we take m ϕ = 5 GeV.
rapidly (see Fig. 5). We find the upper bounds of the mass ratio, beyond which the couplings become non-perturbative (i.e., g µ/τ < √ 4π). Since a larger coupling is required in the cBE approach, a smaller mass ratios are allowed compared to the standard nBE treatment.
FIG. 5: Correct relic density as function of mass ratio and coupling g µ/τ in both nBE and cBE approaches. We take m ϕ = 0.26 GeV for the µ + µ − mode, and m ϕ = 5 GeV for the τ + τ − mode.
Speaking overall, DM relic density in the cBE method is larger than that by using the standard nBE method. When the elastic scattering is strongly suppressed, the temperature of DM particles drops below that of the thermal bath (T χ < T SM ), meaning DM particles do not have enough kinetic energy. The forbidden annihilation thus becomes ineffective due to the inability to overcome the annihilation threshold, leading to a larger abundance in cBE treatment.

IV. PARAMETER SPACE AND EXPERIMENTAL CONSTRAINTS
In this section, we will find the feasible parameter space for the cBE approach in the (m ϕ , g µ,τ ) plane, by requiring the correct DM relic density. The viable parameter space further displays more accurate results compared with nBE. Of course, there are numerous constraints on the model parameters that are imposed by the collider searches and the astrophysical observations across a wide range.
We first study the µ + µ − forbidden channel. The numerical results are shown in Fig. 6. We set g A µ = 0 for the left panel, while g A µ = g µ for the right panel. The different choices of the couplings to muons do not strongly affect the phenomenology.
To elaborate on our findings, we start by searching the parameter boundary of the forbidden annihilation, in which ∆ → 0 and the right relic density are fulfilled. The green line in the plots shows the boundary of the cBE approach; the green shaded region, denoted by the expression ∆(cBE) ⩽ 0, represents the unforbidden space. In the remaining parameter space, where the forbidden annihilations dominate the DM depletion, m χ is chosen at each point to match the correct relic density. The required coupling becomes larger when m µ /m χ > 1, leading to the corresponding curves lying inside the funnel area with a similar shape as the boundary. For comparison, we also repeated the parameter boundary of the standard nBE approach in the same figure, as shown by the gray lines. The gray regions stand for the non-forbidden regions.
The allowed parameter space of the two approaches is noticeably distinctive. The eKD effect reduces the parameter space, compared to that of nBE. Or we can say that to satisfy the relic density requirement, the larger coupling g µ is required for the cBE scenario. The reason is straightforward as already pointed out in the last section. Around the freezeout stage, the temperature of DM decreases resulting in reduced kinetic energy and then inefficient forbidden annihilations. To maintain the DM dilution process, a larger coupling is required for T χ < T .
In the small mediator mass region, the annihilation into pairs of mediators dominates the relic density. We depict it in gray with the label χχ → ϕϕ as the non-forbidden DM region. For experimental constraints, the most crucial parameter is the mediator mass m ϕ . In the following, we show the constraints one by one: • Planck (brown region).
DM annihilations into SM electromagnetically interacting particles could modify the anisotropies of the CMB [47][48][49]. The measurements of the CMB by the Planck satellite [33] can thus put robust constraints on such energy injection processes [35]. In this model, the photon pairs can be produced in DM annihilations via a muon loop, leading to the injection into the CMB. We recast the corresponding constraints from Ref. [23] which are shown in brown.
As pointed out in Ref. [23], when m ϕ near 2m χ which is also close to but smaller than 2m µ , the annihilation cross section of χχ → γγ is enhanced due to σ(χχ . So the CMB constraints exclude the most parameter space for m ϕ < 2m µ and only leave a small part of the allowed room when m ϕ > 2m µ . From Fig. 6, we find that the eKD effects narrow down the allowed parameter space that obtains the right amount of DM relic abundance. Especially the CMB constraints exclude the most parameter space with cBE in the pure scalar interaction scenario (i.e. g A µ = 0).
Secondary muons are produced from the electron beamdump experiments, such as SLAC E137 [36] and Jefferson Lab BDX experiments [40], which can be used to explore the signals of light scalar emission and the muon-scalar coupling via muon-nucleon scattering pro-FIG. 6: Results and constraints on forbidden channel χχ → µ + µ − . The green line represents the parameter boundary of the forbidden DM region with m µ /m χ = 1; and the green shaded region stands for the non-forbidden space that m χ > m µ . The rest space is allowed for the forbidden annihilations to get the correct DM relic density. The gray region labeled ∆(nBE) ⩽ 0 is obtained in the nBE approach. Also, the non-forbidden space that is dominated by the χχ → ϕϕ process is shown. The brown region displays the Planck bounds on the energy injection process [35]. The orange region is excluded by the E137 electron beam-dump experiment [36]. The purple region is excluded by energy loss process of SN1987A [37][38][39]. The projected sensitivities for future beam-dump experiments from BDX [40] (light blue), M 3 [41,42] (dashed blue line), NA62 [43] (red), and NA64-µ [44][45][46] (dashed magenta line) are also shown here.
cess µ + N → µ + N + ϕ. The E137 experiment's null result established exclusion limits on the parameter space, and the upcoming BDX experiment can likewise yield a predicted limit.
In Fig. 6, we display, as shaded areas, constraints from the E137 electron beam-dump experiment [50], and projections from BDX [40]. The constraints exclude part of the available parameter space in the (m ϕ − g µ (g A µ )) plane.
NA62 is a fixed-target experiment at the CERN Super Proton Synchrotron (SPS) that is dedicated to measurements of kaon rare decays, including projected searches on K → µνϕ [43]. Such a decay channel is an excellent probe of new light scalars that couples preferentially to muons. Ref. [51] has derived the probe sensitivity for this process, which can be used to test our parameter space, as shown in red.
Similar to the above beam-dump experiments, NA64µ [44][45][46] and M 3 [41,42] are designed to search light scalars in the muon-nucleon scattering process µN → µNϕ, using muon beams. It's worth noting that the allowed parameter space can be tested in the coming future.
At last, for the parameter space involving new light scalars, we should also consider the constraints from the observation of supernovae cooling. The most famous constraints arise from the energy loss process via the Primakoff effect γp → pϕ in SN1987A [37][38][39]. The excluded region is displayed in purple.
The forbidden annihilation in the di-tau case is discussed here. Annihilation to τ + τ − shares several qualitative features with annihilations to muons. The same computations are performed, including the searches for the forbidden annihilation parameter region, the comparison of cBE and nBE treatments, and the variety of limitations from experimental searches. The results are displayed in Fig. 7 with different parameter settings. Note that the consideration of the experimental constraints is rather different from the χχ → µ + µ − case due to the different mediator mass.
• BaBar (orange region) and Belle II (red region) Searching for e + e − → γ + invisible at the e + e − colliders, such as BaBar and the future Belle II experiments, sets existed and projected constraints on the parameter space. As depicted in Fig. 7, the orange line and region show the constraints arising from BaBar [39,44,52], and the red region shows the exquisite sensitivity from Belle II using 50 ab −1 integrated luminosity [39,54], which indicates that a large portion of the viable parameter space will be tested.
Large Electron Positron Collider (LEP) has set corresponding limits on such channels [44,53] which are shown in brown.
Additionally, there have been several proposals for future Z-factories to search the same processes [55][56][57][58], based on Circular Electron Positron Collider (CEPC) and the Future Circular Collider e + e − (FCC-ee) for instance. The projections of Tera-Z options (with accumulated 10 12 Z's ) provide leading sensitivities in the tens GeV range, shown in purple.
• Planck (blue region) The CMB constraints are similar with the χχ → µ + µ − forbidden channel. The parameter space at m ϕ < 2m τ is much constrained as always, as shown in Fig. 7, due to the enhancement of χχ → γγ when m ϕ ≃ 2m χ .

V. CONCLUSION
In the forbidden DM scenario, consideration of the early kinetic decoupling is not only a correction to the DM relic density but an indispensable ingredient. In this work, we investigate the early kinetic decoupling effect in forbidden channels, where DM annihilation to SM leptons is kinetically forbidden. Specifically we focus on the χχ → µ + µ − and χχ → τ + τ − modes. By analyzing the scattering momentum transfer rate, we found the kinetic equilibrium breaks at about x ≃ 20, which is the same stage of chemical decoupling. So eKD should be taken seriously into account. With different benchmark points, we found that there is a cooling phase during the evolution that causes the DM temperature to deviate from the thermal bath and evolve solely. The decreased kinetic energy of DM particles suppresses the forbidden annihilation rate, which gives rise to larger abundances. The difference between the cBE and traditional nBE methods is significant, within part of the parameter space showing a deviation up to an order of magnitude larger, for both µ + µ − and τ + τ − forbidden channels. We also considered the experimental constraints from beam-dump experiments, collider searches, and astrophysical observations. The viable parameter space in the forbidden DM model has been reduced when using the cBE treatment. Most of the parameter space will be tested by the forthcoming experimental searches.