Enabling Forbidden Dark Matter

The thermal relic density of dark matter is conventionally set by two-body annihilations. We point out that in many simple models, $3 \to 2$ annihilations can play an important role in determining the relic density over a broad range of model parameters. This occurs when the two-body annihilation is kinematically forbidden, but the $3\to 2$ process is allowed; we call this scenario"Not-Forbidden Dark Matter". We illustrate this mechanism for a vector portal dark matter model, showing that for a dark matter mass of $m_\chi \sim \text{MeV - 10 GeV}$, $3 \to 2$ processes not only lead to the observed relic density, but also imply a self-interaction cross section that can solve the cusp/core problem. This can be accomplished while remaining consistent with stringent CMB constraints on light dark matter, and can potentially be discovered at future direct detection experiments.


I. INTRODUCTION
The particle physics nature of dark matter (DM) is still a mystery despite undeniable evidence of its gravitational interactions. The observed relic abundance of DM may provide a clue to its non-gravitational interactions, as in the classic weakly interacting massive particle scenario, where the freezeout of 2 → 2 annihilation of DM particles to the standard model (SM) particles sets the late-time abundance of DM. Many variations on the standard thermal freezeout scenario have recently been considered (e.g. ); in this article, we point out that even for simple and weakly-coupled dark sectors, 3 → 2 annihilationsas illustrated in fig. 1 -can play a critical role.
For weakly-coupled DM, 3 → 2 processes are usually considered to be subdominant to their 2 → 2 counterparts at the time of freezeout, but if the latter are kinematically suppressed while 3 → 2 is unsuppressed, the situation is more complex. This can occur when the DM couples to a "mediator" particle with a mass somewhat larger than that of the DM itself, as might arise in a hidden sector characterized by a single scale.
Kinematic suppression of 2 → 2 annihilation, leading to a novel cosmological history during freezeout, was previously invoked in the "Forbidden DM" [9] and "Impeded DM" [25] scenarios; the new feature in our study is the presence of a kinematically allowed dark-sector 3 → 2 annihilation channel. We refer to this scenario as Not-Forbidden Dark Matter (NFDM). The 3 → 2 channel is also important in the Strongly Interacting Massive Particle (SIMP) scenario [3], but work on SIMPs has focused on strongly coupled theories with scalar DM [4,26], whereas NFDM is a more generic mechanism: it is potentially important in any situation where 2 → 2 annihilations within the dark sector are kinematically suppressed, and has no obvious dependence on whether the * jcline@physics.mcgill.ca † hongwan@mit.edu ‡ tslatyer@mit.edu § weixue@mit.edu DM is fermionic or bosonic or whether the dark sector coupling is strong or weak. Hidden sector or multicomponent DM models may have regions of parameter space where NFDM is an important mechanism to consider. We illustrate our paradigm with a Dirac fermion DM charged under a hidden U (1) symmetry, with dark gauge boson A . This mediator can provide a portal to the SM by having a small coupling to the electromagnetic current J µ EM through a kinetic mixing term ( /2)F µν F µν . In the mass basis, the Lagrangian becomes The gauge coupling is α = g 2 /4π, and / D ≡ / ∂ − ig / A . It is clear in this basis that there is no tree-level coupling between χ and the SM photon. We can also consistently assume that the dark Higgs boson giving mass to A is arXiv:1702.07716v2 [hep-ph] 13 Oct 2017 very heavy and can be neglected in the effective description [27]. Depending upon the size of the kinetic mixing parameter , there are two possible regimes of interest: (1) is relatively large, such that the hidden sector and the SM sector have the same temperature before DM freezeout, while is still small enough so that 3 → 2 and 2 → 2 reactions involving only hidden sector particles dominate over annihilation of χ to SM particles; (2) For sufficiently small 10 −8 , the hidden sector will have its own temperature and in the limit → 0, it becomes secluded: both χ and A contribute to the ultimate DM density.
In section II, we discuss the freezeout history of the NFDM model, and by solving the Boltzmann equations, we determine the dark sector parameter values {m χ , m A , } that are consistent with the observed relic density. In section III we incorporate constraints from a variety of astrophysical and laboratory searches, showing that a significant parameter region is allowed while realizing the NFDM mechanism. Conclusions are given in section IV. In the appendix, we present a more detailed account of how the order of freezeout of the various reactions determines the relic abundance; the dependence of our results on the temperature of the dark sector; how the constraints change with m A /m χ , and cross sections for the relevant scattering processes.

II. COSMOLOGY
Previous studies of the vector-portal DM model, shown in eq. (1), have divided the parameter space into two broad regions: m χ < m A or m χ > m A . In the latter case, the dominant process at the epoch of thermal freezeout is χχ → A A followed by A decays to SM particles, whereas when m χ < m A , the s-channel annihilation χχ → ff to SM particles f via off-shell A is dominant. This regime is ruled out for m χ ∼ MeV-GeV by CMB constraints [28][29][30][31][32].
In the present work, however, we are interested in the intermediate region m χ m A , where it is possible for the 3 → 2 scatterings χχχ → χA or χχA → χχ to have an important effect on the dark matter abundance. The system is governed by the coupled Boltzmann equations for the χ and A densities. For m χ m A , the relevant terms in these equations are where n χ (n χ,0 ) denotes the (equilibrium) density of χ +χ, and similarly n A (n A ,0 ) for the dark photon. Throughout this paper, we have assumed zero chemical potential for χ and χ, and take the densities of χ and χ to be equal. The 1/4 in the first term of eq. (2) is the symmetry factor for Dirac DM, taking into account the two identical particles in the initial state and the fact that each annihilation process removes a χχ pair. The conjugate process χχχ → χA is also accounted for in this factor. The relative numerical factors between the two equations are consistent with the way each process changes the number density of χ and A ; for example, the factor of 1/4 and 1/8 in the first terms of eq. (2) and (3) respectively are consistent with the fact that the 3 → 2 process has a net effect of removing a χχ pair and producing a single A . A detailed discussion of the derivation of the Boltzmann equation for 3 → 2 processes can be found in [33].
Other 3 → 2 processes such as χχA → χχ, 3A → χχ etc. are important only in the case of m A /m χ < 1 and = 0 in Sec. II B. The complete Boltzmann equations containing all of these processes are shown in eq. (D1) and (D2) in Appendix D. All numerical results in this paper across the full range of m A /m χ considered were obtained using the complete equations. Expressions for the cross sections are given in Appendix E.
We will focus on the two regimes where (1) the hidden sector and the SM remain in thermal equilibrium, requiring values of the kinetic mixing 10 −7 (but still small enough to avoid dominance of the χχ → e + e − process); (2) the hidden sector is secluded from the SM, → 0.
A. Kinetic equilibrium with the SM For sufficiently large , the scattering process χe ± → χe ± is fast enough to keep the dark and visible sectors in kinetic equilibrium, T d = T SM . By comparing the rate inferred from the χe ± → χe ± cross section to the Hubble rate H at DM freezeout, we estimate the condition to be taking x f ∼ 20, the e ± to be relativistic, and m A m χ . This leaves a significant range of 10 −8 − 10 −4 , depending upon m χ , such that A -mediated annihilations χχ → e + e − are out of equilibrium (a requirement of our scenario), as will be shown below. We take the dark sector masses to be in the ranges m χ ∼ < 10 GeV and m χ ∼ < m A < 2m χ . The lower bound on m A makes χχ → A A kinematically inaccessible, while the upper bound forbids the A → χχ decay channel. If m A > 2m χ , the number-changing process χχχ → A χ effectively becomes number-conserving, χχχ → χχχ. In terms of the parameter r ≡ m A /m χ , the relevant range is thus 1 r 2.
The Boltzmann equations are solved numerically, and the results shown in fig 2. As an example, fig. 2(a) illustrates the evolution of the χ and A abundances as a function of x ≡ m χ /T with m χ = 0.2 GeV, gauge coupling α = 1, kinetic mixing = 10 −6 and the ratio r = 1.9. This example has been chosen to emphasize the importance of 3 → 2 scatterings, but similar results are obtained for r 1.5. Here, in the case with only 2 → 2 annihilation, the DM abundance would reach its relic value at x f ∼ 20; in our NFDM case, in contrast, the 3 → 2 processes and decay of the A control the freezeout, and their interplay leads to an extended freezeout con-tinuing to x f ∼ 60. If we neglect the 3 → 2 process the resulting abundance is overestimated by several orders of magnitude. It is noteworthy that Y A departs from the equilibrium abundance at late times, even though the rate for A → e + e − exceeds the Hubble rate, because the 3 → 2 or 2 → 2 processes can also strongly affect the evolution of n A .
In fig. 2(b) we plot the contours in the m χ -r plane matching the observed relic density [31], for several values of α and . We consider values of α ≤ 4π, since every loop integral introduced in a Feynman diagram typically introduces an additional factor of α/4π, and so perturbativity is naively maintained for this range of α . n A = n A ,0 correponds to large , where the rate for A → e + e − dominates the rates for either of the two annihilation processes that generate A s. The region r 1.5 corresponds to the Forbidden DM regime, and ref. [9] studied this regime with the assumption of n A = n A ,0 : smaller values of show increasing deviation from the relic density contours obtained from this assumption, even for r < 1.5. For the rest of the paper, we will focus on the NFDM region 1.5 r < 2, where the 3 → 2 process leads to a strong transition in the behavior of the relic density contour, with the exact value of r for the transition depending on the coupling constant α .
Normally the DM relic density is set by the strongest annihilation channel, which is also the last to freeze out, since only a single Boltzmann equation for DM is considered. This applies when is large, forcing n A n A ,0 (dashed contours). These contours turn to the right as r → 2 because the 3 → 2 cross section diverges, σv 2 χχχ Thus obtaining the correct relic density as r → 2 requires a larger value of m χ .
In contrast, for moderate values of , the NFDM mechanism applies, where the two coupled Boltzmann equations must be solved together. In general, we find that typically the two strongest processes (either annihilations or decays) keep the coupled system in equilibrium until the rate for one process (per χ particle) becomes comparable to the Hubble rate, and thus any weaker processes are not relevant for determining the relic abundance. In this regime, typically the decay of A → e + e − and either the 3 → 2 or 2 → 2 annihilation are the relevant processes. In particular, for r 1.5 − 1.8, the 3 → 2 scatterings are faster than 2 → 2, and so they dominate the freezeout, as shown in fig. 2(a). The combination of 3 → 2 scatterings and A decays can lead to a non-equilibrium density for the A particles during the freezeout of the 3 → 2 process if is sufficiently small (e.g. ∼ 10 −6 − 10 −7 ), resulting in a lengthy freezeout and an -dependent relic density. This behavior corresponds to the divergence of the dashed and solid contours in fig. 2(b) at large r.

B. Secluded hidden sector
Next we consider the limit of → 0, so that the dark photon is effectively stable, and the hidden sector is secluded. This analysis can be easily applied to multicomponent DM models. Even though secluded hidden sectors are in general difficult to probe due to the lack of any interaction with the SM, they are not entirely impossible to study. Secluded hidden sectors can, for example, be constrained by the number of relativistic degrees of freedom during Big Bang Nucleosynthesis (BBN). Furthermore, in the U(1) theory considered here, the relic abundance is set by the coupling strength α , which in turn determines the self-interaction cross section in the dark sector. This cross section is a prediction of the model, and has observable consequences for structure formation, which can in principle be highly constraining.
Moreover, the secluded case is a useful limit that gives insight into the region of parameter space where is small but non-zero, so that kinetic equilibrium cannot be maintained with the SM. Despite the small couplings to the SM, this regime can still be effectively probed by observations of the cooling of SN1987a [34]. The secluded limit is also highly instructive as an illustration of the rich behavior that can occur in the U (1) vector portal DM model when the 2 → 2 and 3 → 2 annihilations are the dominant processes at freezeout.
To avoid warm or hot dark matter [12], we assume that χ couples additionally to some relativistic degree of freedom φ until freezeout, strongly enough to maintain thermal equilibrium in the dark sector so that the DM temperature redshifts with the Hubble expansion in the conventional manner, T ∼ 1/a. However, the coupling of φ to χ should be sufficiently weak that annihilation of χχ → φφ is negligible during freezeout, to make the NFDM freezeout mechanism dominate over conventional 2 → 2 annihilation. For a concrete model of how this can be achieved, we take φ to be a light scalar charged under some additional U (1) symmetry, interacting with the dark sector through the dimension-5 operator (1/Λ)χχφ * φ. The T ∼ 1/a dependence is maintained by χφ → χφ scatterings, which has a rate that scales as n φ σv χφ→χφ , while the χχ → φ * φ rate scales as n χ σv χχ→φ * φ . To obtain a parametric estimate for a value of Λ that would maintain both T ∼ 1/a and subdominance to the 2 → 2 and 3 → 2 processes considered in eq. (2) and (3), we take σv χφ→χφ ∼ σv χχ→φ * φ ∼ 1/Λ 2 , and look for values of Λ which ensure that the χχ → φ * φ annihilation rate is subdominant up to the point of freeze-out of the two main processes. This condition is most difficult to satisfy in the case where r = 2, and the 2 → 2 rate becomes highly suppressed. Nevertheless, we find that in this limit, a suitable range of Λ is m m χ m pl , which for GeV dark matter corresponds to 10 6 Λ/GeV 10 9 .
More generally, the dark sector has its own temperature T d which need not be the same as that of the visible sector, T SM ; it is determined by details of the thermal cosmological history such as the efficiency of reheating into the dark sector after inflation. The relic abundance in this case depends upon the unknown parameter γ ≡ T d /T SM , but in a simple way: Y χ ∝ γ p(r) where p(r) ∼ 1.6 − 1.8 depends upon the mass ratio r = m A /m χ . Here we illustrate the case of γ = 1.
The evolution of n χ and n A for the secluded dark sector is shown in fig. 3(a), taking m χ = 35 MeV, α = 1, r = 1.95 as an example to illustrate the important interplay between the 2 → 2 and 3 → 2 interactions. Keeping only the 3 → 2 reaction would predict that A becomes the dominant DM component, whereas in reality it remains highly subdominant. Again the freezeout process is prolonged, starting with the decoupling of 2 → 2 scatterings at x ∼ 20, while the 3 → 2 reactions decouple at x ∼ 150. Interestingly, n A temporarily grows between these two times, allowing the 2 → 2 rate to come back above H just before freezeout completes.
In fig. 3(b), we plot contours corresponding to the observed thermal relic density in the m χ -r plane for different values of α . In the following we give a brief explanation of the contour shapes in the regions r 1, 1 r 1.5 and 1.5 r 2, which each show a distinct qualitative behavior: (1) r 1. Being lighter than χ, A is the dominant DM constituent. The fastest process in this mass range is the 2 → 2 process χχ → A A . Significantly below r = 1, the second fastest process is 3A → χχ, since n A ,0 > n χ,0 . Near the threshold, with n A ,0 ∼ n χ,0 , all of the other possible 3 → 2 processes (χχA → χχ, χχA → χχ, 10 100 10 -15 10 -10 χχA → A A , χA A → χA , as well as χχχ → χA plus any conjugate processes) become important. The relic abundance curves in fig. 3 are computed with all of these processes taken into account in the complete Boltzmann equations shown in eq. (D1) and (D2).
(2) 1 ∼ < r ∼ < 1.5. χ is the dominant DM component. The fastest reaction is A A ↔ χχ, and it enforces n A = n A ,0 n χ /n χ,0 during the freezeout, and the second fastest reaction is now χχχ → χA , which determines the DM abundance. The 3 → 2 rate goes as n 2 χ σv 2 , which depends only weakly on r through the phase space. Therefore there is no strong correlation between the abundance and r in this region.
(3) 1.5 ∼ < r ∼ < 2. χ is the dominant DM component, but now its abundance is determined by the two freezeout events A A → χχ (whose rate becomes comparable to Hubble at later times) followed by χχχ → χA . At large r 2, just before freezeout completes, both reactions are faster than H, allowing one to estimate the freezeout times. Taking the 2 → 2 and 3 → 2 rates ∼ H, and n A n A ,0 n χ /n χ,0 enforced by fast 3 → 2 scatterings, we can analytically derive contours consistent with the numerical results.

III. CONSTRAINTS
The parameter space of NFDM is constrained by a variety of experimental observations: i) dark photon limits coming from the cooling of SN1987a [34]; ii) similar bounds from beam dump experiments [35,36]; iii) limits on the thermally-averaged cross section of χχ → e + e − deduced from the CMB power spectrum measured by Planck [31,32,[43][44][45], and iv) direct detection constraints on the dark matter-nucleon scattering cross section from PandaX-II [37], LUX [38] and CDMSLite [39]. Although we have only assumed a coupling to electrons in much of this analysis for simplicity, these direct detection limits are relevant to the vector-portal DM model considered here.
Future direct detection experiments including Super-CDMS SNOLAB [40], as well as electron scattering off germanium [42,[46][47][48] and graphene [41] are also shown in the same plot. Other current limits from XENON10 [49], indirect detection [50] and lower bounds on m χ from N eff [51] are sub-dominant to the current constraints presented here and are not shown. Fig. 4 summarizes these constraints in the m χ -plane for the illustrative value of r = 1.8, with α fixed to give the correct present-day relic density, subject to the perturbativity constraint α ≤ 4π. At a large and small m χ , conventional freezeout from χχ → e + e − annihilations dominates over the NFDM mechanism, but this is ruled out by the CMB constraint. The approximately horizontal red dashed contour shows the minimum value of for which the visible and dark sectors are in kinetic equilibrium, estimated in eq. (4).
Self-interactions between dark matter particles with a cross section σ SI ∼ 0.1 σ SI /m χ 1 cm 2 g −1 can potentially resolve the core-cusp and the too-big-to-fail  [35,36] (pale orange), SN1987a cooling [34] (blue), direct detection [37][38][39] (yellow) and perturbativity, α ≥ 4π (gray). The projected reach of SuperCDMS [40] (orange dot-dashed line), electron ionization of graphene [41] (magenta dot-dashed line) and germanium in a low-threshold experiment [42] (blue dot-dashed line) are also shown. The curve near ∼ 10 −7 indicates where kinetic equilibrium with SM is established (red dashed line). The region of parameter space where the self-interaction cross section exceeds current limits (σ/mχ > 1 cm 2 /g) (purple), and the region where the self-interaction cross section can potentially solve the smallscale structure problems (0.1 cm 2 /g < σ/mχ < 1 cm 2 /g) (purple dashed lines) are displayed. The purple arrow points into the region allowed by self-interaction bounds, above and to the right of the line. The A decay rate is faster than H at freezeout above the lowest (blue) curve.
problems of small-structure formation with cold DM [52][53][54] while remaining consistent with experimental constraints, which set an upper bound of between 1 -2 cm 2 g −1 [55][56][57]. A DM mass of m χ ∼ (0.1 − 1) GeV with ∼ 10 −7 − 10 −6 in our model leads to a velocityindependent self-interaction cross section that lies within this range, and can provide a possible solution to both puzzles (though recent analysis of clusters indicates some preference for a velocity-dependent cross section [58]). The preferred region is between the purple dashed lines in fig. 4, while the cosmologically constrained region is shown in purple.

IV. SUMMARY AND OUTLOOK
We have demonstrated a novel scenario called Not-Forbidden Dark Matter, where an allowed 3 → 2 annihilation process compensates for its conventional 2 → 2 counterpart being kinematically forbidden during thermal freezeout. This mechanism can be potentially important in a variety of hidden sector models, including vector-portal, scalar-portal and composite DM. The DM mass and the mediator (or second DM) mass are of the same order, which would naturally arise in a hidden sector characterized by a single scale.
Taking the vector-portal DM model as an example, we found that in some parts of the NFDM parameter space, the combined effect of 3 → 2, 2 → 2 and A decay channels is to significantly prolong the period of freezeout. The commonly-neglected 3 → 2 annihilation channel can change the predicted relic density by orders of magnitude. Although this model is restricted by an abundance of experimental constraints, viable examples remain in the mass range ∼ (0.1 − 1) GeV, with a self-interaction cross section that is coincidentally of the right order for solving the small scale structure problems of ΛCDM cosmological simulations. This is a well-motivated target for future direct detection [40,46,47] and dark photon searches [59][60][61][62][63].
While we were completing this work, [64] appeared, presenting a related idea. Their work focuses on keV-MeV scalar DM and requires additional scalar "assister" particles.
As mentioned above, an essential difference between NFDM and conventional DM freezeout is the importance of tracking the evolution of both the DM χ and the mediator particle (in our model, A ), by solving the coupled Boltzmann equations (eq. (2) and (3) for relevant terms when r 1, eq. (D1) and (D2) for the complete equations) for their respective densities. The presence of two equations implies that more than one scattering (or decay) process can be important for determining the final abundance; hence both the fastest and second fastest reactions are typically relevant. This is in contrast to conventional DM freezeout based upon a single Boltzmann equation, where the abundance depends upon the strongest channel. In the large limit of our model, A decay is the fastest process, and enforces equilibrium of A , n A = n A ,0 . Hence smaller values of are necessary to realize the rich cosmology that comes from the interplay of the coupled Boltzmann equations of χ and A . To simplify the subsequent discussion, we assume that these -suppressed reactions are negligibly slow, i.e. we work in the secluded dark sector regime of the NFDM model.
It is useful to define the net rate of 3 ↔ 2 or 2 ↔ 2 interactions per χ or A particle by considering the collision terms in the Boltzmann equations, written in the Likewise, we define 2(n A /n χ )R A (χχχ → χA ) ≡ R χ (χχχ → χA ) and so on for the unidirectional rates. In this way, the signs for these definitions have been chosen so that all of the rates of individual sub-processes are now positive, although the overall rates R χ and R A can have any sign. When m A > m χ and T < m χ , m A , the lower density of A relative to χ implies that R χ (3 ↔ 2) (R χ (2 ↔ 2)) is generally smaller in magnitude than R A (3 ↔ 2) (R A (2 ↔ 2)). Thus the rates R χ tend to fall below H earlier than the corresponding rates R A . This separation between freezeout of χ and A is the origin of the prolonged duration of the overall freezeout process. Suppose that only one channel, for example 2 → 2, occurs fast enough such that R χ (A A →χχ) H; then this rate tends to be nearly equal to that of the reverse reaction, R χ (χχ → A A ), enforcing the condition n 2 A n 2 A ,0 n 2 χ /n 2 χ,0 (though the cancellation is imperfect, so that the total rate R χ (2 ↔ 2) is also typically greater than H). This by itself is not sufficient to force both the χ and A densities to track their equilibrium values. For that, one generically needs both R χ (3 ↔ 2) > H and R χ (2 ↔ 2) > H so that both independent combinations n χ −n χ,0 and n A −n A ,0 are driven to zero. 1 This is 1 The typical behavior is that the strongest process is such that always true at sufficiently early times, allowing us to use equilibrium initial conditions for the Boltzmann equations. The DM density n χ starts to deviate from equilibrium when the rate of the weaker annihilation channel becomes comparable to H; hence the second-strongest channel initiates the freezeout process.
To illustrate this behavior, we show some examples of the evolution of the rates in fig. 5a, 5b and 6. Each example has the same DM mass m χ = 70 MeV, coupling α = 1, and kinetic mixing = 0, but different values of r = 1.4, 1.7, 1.9. In these figures, the dot-dashed lines corresponding to the evolution of DM number density are shown to highlight the time of DM freezeout. For r = 1.4, the freezeout period is relatively short; for r = 1.7, freezeout is prolonged; and for r = 1.9, the freezeout is prolonged further and may indeed be thought of as two separated freezeouts.
In fig. 5b, we show the two rates R A (3 ↔ 2) and R A (2 ↔ 2), which are much larger than H; these cancel each other to order H. The behavior is similar for other values of r.
. This relation is demonstrated in fig. 5a and 6.
Comparing these net rates however does not tell us which process controls freezeout. Instead, we should look at the dashed lines, which indicate the unidirectional rates, R χ (A A →χχ) and R χ (χχχ → χA ). Processes are out of equilibrium when these dashed lines overlap with the solid lines. From the unidirectional rates, we can identify the weaker annihilation channel and thus which process initiates the freezeout. For r = 1.4, the weaker process is 3 → 2, and fig. 5a confirms that the freezeout is indeed triggered by 3 → 2. For r = 1.7 and r = 1.9, the dashed line for R χ (A A →χχ) in fig. 6 merges with the solid line, R χ (2 ↔ 2), when the rate is about 3H. It is the weaker channel 2 → 2 that initiates freezeout.
One difference between r < 1.5 and r > 1.5 in figs. 5a versus 6 is that r > 1.5 normally has a longer freezeout. The duration depends upon whether the rate of the weaker annihilation channel is sensitive to n A . For r > 1.5, the weaker process 2 → 2 has the rate R χ (A A → χχ) ∼ σv A A →χχ n 2 A /n χ . Prior to the final freezeout, the larger 3 → 2 rate imposes the constraint that n A n A ,0 n 2 χ /n 2 χ,0 ∼ r 3/2 x 3/2 m −3 exp((2 − r)x)n 2 χ . Since n A increases with time, this means that the 2 → 2 rate R χ (A A →χχ) can be kept at the same order as H for a long period. For this reason, the freezeout is prolonged. For r < 1.5, the duration is relatively short, because the rate of the weaker 3 → 2 process goes as R χ (χχχ → χA ) ∝ n 2 χ , where n χ is decreasing with time. Armed with our insight that the second-strongest channel matters critically for freezeout, and having unboth the forward and backward rates exceed H, as well as their difference. For the second-strongest, only one of these need be greater than H.
derstood the reason that the freezeout process is longer for r 1.5, we can also explain the shape of the contours in fig. 2(b) in the main text. As discussed briefly in the main text, there are several important regimes: 1. For n A = n A ,0 , corresponding to large , the two Boltzmann equations are reduced to one, and the shape of the contours can be understood using the usual parametrics of DM freezeout. This behavior occurs for the contours overlapping the dashed contours in fig. 2(b).
2. For = 0, we recover the secluded case discussed above, where the interplay of the 3 ↔ 2 and 2 ↔ 2 processes controls the freezeout. This behavior also occurs in the region where r ∼ 1.5 and α = 10, because the 3 → 2 and 2 → 2 rates for DM are significantly larger than R χ (A → e + e − ) ≡ Γ(A → e + e − )n A /n χ , so that Γ A →e + e − can be neglected. 2 3. For moderate , the rates of the three processes should be compared in order to determine which is weakest, and hence irrelevant to the DM freezeout.
The relevant rates to compare are R χ (χχχ → A χ), R χ (A A → χχ) and R χ (A → e + e − ), evaluated at the Hubble crossing time of the annihilation processes. Consider the case where 2 → 2 has a lower rate than 3 → 2, such that it falls below H first.
Whether the DM density freezes out or not at this time depends on the relative sizes of R χ (A → e + e − ) and R χ (A A → χχ). When R χ (A A → χχ) < R χ (A → e + e − ), the larger R χ (A → e + e − ) rate in the coupled Boltzmann equations provides enough constraints to keep n χ and n A near their equilibrium values. An example of this more complex case is shown in fig. 2(a) of the main text, where the freezeout starts when R χ (3 ↔ 2) ∼ H. Using the Boltzmann equation of A , this rate is determined by More broadly, this case is realized when = 10 −6 , and either r is close to 2, or α is large and r > 1.5. Fig. 2(b) shows that in this region the = 10 −6 contours (solid lines) do not overlap with the dashed contours, for which the constraint n A = n A ,0 is imposed. In this regime the 3 → 2 rate is the largest, and when R χ (A A → χχ) ∼ H, R χ (A → e + e − ) > R χ (A A → χχ). The freezeout is thus controlled by 3 → 2 processes and the decay of A . In this case the A decay rate is not fast enough to keep the A abundance in equilibrium, and both n A and n χ are increased during freezeout relative to their values when the A s remain in equilibrium. The χ annihilation rate thus needs to be increased to maintain the correct relic density, requiring lower χ masses; this is the reason that the contours in fig. 2(b) bend toward lower masses as is decreased, for large r. and ultimately the DM relic abundance. In the Boltzmann equations (eq. (2, 3) in the main text), taking T d = T SM changes the equilibrium densities, so that n χ,0 ∼ exp(−m χ /T d ) = exp(−x/γ), where we have defined γ ≡ T d /T SM , and x is still given by x ≡ m χ /T SM . Likewise, n A ,0 ∼ exp(−rx/γ). Keeping in mind that H is determined by T SM , we can solve the Boltzmann equations for n χ (x) and n A (x) with the γ-dependence coming from the equilibrium densities. Fig. 7 shows the behavior of the ratio of relic abundances Ω c (T d )/Ω c (T d = T SM ) as a function of T d for 0.1 ≤ γ ≤ 1. Having T d < T SM leads to an earlier freezeout, since the exponential decrease in the equilibrium densities occurs more rapidly. For values of r where the backward and forward 3 → 2 processes fall out of equilibrium at freezeout, we expect that n 2 χ ∼ H/ σv 2 χχχ→χA ∼ 1/x 2 f , and therefore that the relic abundance scales as Ω c ∼ x 2 f . For 1 < r 1.5 where the 3 → 2 process determines the DM abundance, the exponential dependence of n χ,0 with x/γ results in Ω c ∼ γ 2 . On the other hand, in the case of r 2, the second freezeout occurs at n χ ∝ n 4 χ,0 /n 2 A ,0 , and a similar argument leads again to Ω c ∼ γ 2 . At intermediate values of r, both the 3 → 2 and 2 → 2 processes freeze out at similar times. For a 2 → 2 freezeout, n χ ∼ H/ σv χχ→A A , and as a result Ω c ∼ γ. Qualitatively, we expect the γ dependence to lie between these two regimes for intermediate values of r.   fig. 4 in the main text, but also exhibit several distinct characteristics that we explain here.
At r = 1.4, the transition from the secluded limit ( → 0) to the kinetic equilibrium limit occurs in the range ∼ 10 −9 − 10 −6 , which leads to a rapid decrease in α between these two phases at a fixed value of m χ . This explains the rapid change in the behavior of the region with suitable self-interaction for m χ 100 MeV.
At r = 1.6, the most distinctive feature is the change in behavior of the region where annihilations to e + e − dominates at m χ ∼ 100 MeV. At masses smaller than this point, the 2 → 2 process is freezes out last, while at larger masses, it is the 3 → 2 process which does so. This difference accounts for the change in the slope of the boundary. There is no such transition for the other two cases, since at r = 1.8, the 3 → 2 process always freezes out last, while for r = 1.4 it is the 2 → 2 process instead.
For all values of r, a significant part of the m χ − parameter space is still consistent with the present-day relic density while evading experimental constraints, showing that the NFDM scenario is robust against taking differ-ent values of r 1.5.

Appendix D: Complete Boltzmann equations
The complete Boltzmann equations, including all relevant 2 → 2 and 3 → 2 processes for the full range of r considered is: The symmetry factors preceding each term properly account for the number of identical particles in the initial state, the net number of particles created or destroyed in each annihilation process, as well as conjugate processes. These equations are used in all numerical calculations shown in the paper.

Appendix E: Cross sections and decay rates
The decay rate for A → e + e − is For scattering cross sections, the thermally-averaged 2 → 2 cross section for the process 1 + 2 → 3 + 4 is given by where g i is the number of degrees of freedom and f i is the phase space distribution of species i. The averaged squared matrix element |M| 2 is averaged over both the initial and final state degrees of freedom. S f = i n i ! is a symmetry factor, where n i is the number of identical particles of species i in the final state. Initial state symmetry factors are included explicitly in the Boltzmann equation, eq. (2) and (3). This convention may differ from other sources in the literature.
Similarly, the thermally-averaged 3 → 2 cross section for the process 1 + 2 + 3 → 4 + 5 is where λ(x, y, z) ≡ (1−(z +y) 2 /x 2 )(1−(z −y) 2 /x 2 ). This expression agrees with the result for the specific process of 3χ → 2χ computed in [67]. 3 In Table I, we list all of the number changing processes that are included in the Boltzmann equations eq. (D1) and (D2), the initial-and final-state averaged squared matrix element |M| 2 of each process as well as the phase space factor P such that σv or σv 2 = P |M| 2 . We define r ≡ m A /m χ throughout.
Two other processes that are important to our analysis are χe ± → χe ± which maintains kinetic equilibrium between the dark sector and the SM, and dark matter-dark matter scattering.