{'fraction_non_alphanumeric': 0.03935306187730645, 'fraction_numerical': 0.016992658360185063, 'mean_word_length': 4.295789778126679, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Enhanced n-body annihilation of dark matter and its indirect signatures

We examine the possible indirect signatures of dark matter annihilation processes with a non-standard scaling with the dark matter density, and in particular the case where more than two dark matter particles participate in the annihilation process. We point out that such processes can be strongly enhanced at low velocities without violating unitarity, similar to Sommerfeld enhancement in the standard case of two-body annihilation, potentially leading to visible signals in indirect searches. We study in detail the impact of such multi-body annihilations on the ionization history of the universe and consequently the cosmic microwave background, and find that unlike in the two-body case, the dominant signal can naturally arise from the end of the cosmic dark ages, after the onset of structure formation. We examine the complementary constraints from the Galactic Center, Galactic halo, and galaxy clusters, and outline the circumstances under which each search would give rise to the strongest constraints. We also show that if there is a population of ultra-compact dense dark matter clumps present in the Milky Way with sufficiently steep density profile, then it might be possible to detect point sources illuminated by multi-body annihilation, even if there is no large low-velocity enhancement. Finally, we provide a case study of a model where 3-body annihilation dominates the freezeout process, and in particular the resonant regime where a large low-velocity enhancement is naturally generated.


Introduction
In recent years there has been great interest in dark matter (DM) models with modified thermal histories, where DM annihilation is still responsible for setting the late-time abundance of DM, but the annihilation involves three or more DM (or DM-like) particles (e.g. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]). This pushes the preferred mass scale for the DM below 1 GeV [2], a range of interest for novel direct detection searches (see e.g. [18] for a recent review) and which has interesting implications for DM self-interaction rates, while simultaneously evading indirect-detection constraints on sub-GeV DM whose abundance is set by ordinary two-body annihilation [19]. Such light DM models are quite challenging to probe with conventional direct detection and at high-energy colliders. Indirect detection signals could be present if the final state of the annihilation process involves Standard Model (SM) particles, or unstable particles which subsequently decay back to the SM; this feature is not present in all models with higher-body annihilation processes, but even when it is, the conventional wisdom is that these signals will always be unobservably tiny long after freezeout, since they are suppressed by an extra power of the DM density compared to the usual case of 2-body annihilation.
There are several possible caveats to this argument: • For parameter regions and search channels where 2-body thermal relic cross sections are already very strongly excluded, it might be possible to observe a 3-body thermal relic signal.
• In ultracompact minihalos or DM cusps around black holes, where the DM density could potentially be very large, the suppression of n-body processes relative to 2-body would be reduced, and it might be possible to observe a signal.
• The DM annihilation rate may depend on parameters other than the DM density; in particular, the cross section may be a function of the relative velocity of the interacting particles, and could be strongly enhanced at low velocities, with larger enhancements for n-body (for n > 2) annihilation compared to 2-body annihilation.
In this work, we will explore these three caveats, to address the question of whether it could ever be possible to detect SM particles produced by n-body annihilation in such models, with n > 2, in cases where (at least some of) the immediate products of annihilation decay back to the SM if they are not SM particles themselves.
In particular, we note that the partial-wave unitarity bound on the annihilation rate has a much stronger velocity dependence as n increases, permitting annihilation rates scaling as rapidly as v −(1+3(n−2)) . In the presence of effective long-range interactions, multiple partial waves can contribute and the velocity scaling can be even stronger. We present an example model of a secluded dark sector where due to the presence of a resonance, the annihilation rate has a very strong velocity scaling down to some regulating velocity scale, leading to greatly enhanced signatures in low-redshift indirect probes (resonant 3 → 2 interactions were also considered e.g. in Ref. [9]). The phenomenology of such models is similar in some ways to Sommerfeld-enhanced [20] or Breit-Wigner-enhanced [21] DM annihilation in the two-body case, but the different density scaling and the much sharper velocity dependence necessitates a reconsideration of indirect-detection limits.
We will explore indirect-detection limits on DM at the keV-GeV scale from the cosmic microwave background and X-ray/γ-ray observations of galaxies and clusters. We will attempt to proceed in a maximally model-independent manner, estimating limits on the parameter which describes the amount of power injected through such annihilation processes. This will allow us to compare the constraining power of different search channels in a largely modelindependent way. We focus here on indirect detection as it constrains the n-body annihilation process directly.
Models of this type could also potentially be detectable at colliders or in direct detection experiments, depending on the strength of the coupling between the dark sector and the SM; indirect detection signals are relatively insensitive to this coupling if the dominant annihilation channel involves only dark sector particles, as the coupling only controls the decay lifetime of unstable dark-sector particles into the SM (which can be quite long without affecting the constraints). However, appreciable couplings between the dark sector and the Standard Model will tend to generate observable signals from 2-body annihilation (directly into SM particles) that overwhelm any n-body signal for n > 2, so we are driven to consider models with very small SM couplings that in turn have very suppressed signatures in terrestrial experiments.
The paper is organized as follows: in Section 2, we provide general parametric estimates for freezeout through velocity-independent n-body processes, the strength of such n-body processes relative to 2-body annihilation in regimes relevant for indirect detection, and upper bounds on the n-body annihilation signal arising from upper bounds on the DM density. In Section 3, we discuss upper bounds from unitarity on the n-body annihilation cross section at low velocities, for general n, and outline some implications if a large low-velocity enhancement is present. In Section 4, we discuss our methodology for estimating contributions to annihilation signals from DM halos and subhalos; this is already a major source of uncertainty for 2-body annihilation and we will see that the effect is amplified for n-body annihilations with n > 2. In Section 5, we derive and compare the sensitivity of various indirect searches to 3and 4-body annihilation, attempting to work as model-independently as possible. In Section 6, we discuss the potential for detecting annihilation signals from a population of UCMHs. In Section 7, we perform a case study of one specific class of models, where 3-body annihilation dominates freezeout and there is potential for large low-velocity enhancements. We present our conclusions in Section 8.
2 Parametric estimates for velocity-independent n-body annihilation

Cross sections for n-body processes
Let us first briefly discuss notation. We will use the notation σv n−1 to denote the annihilation rate coefficient for a given velocity distribution of n annihilating particles. If the n annihilating particles are all distinguishable, and have mass densities ρ 1 , · · · , ρ n , and masses m χ 1 , · · · , m χn (here we assume all initial-state particles are non-relativistic), then the annihilation rate per unit volume per unit time is σv n−1 (ρ 1 /m χ 1 ) · · · (ρ n /m χn ). If some of the annihilating particles are identical, this expression must be corrected by a combinatoric factor; each set of j identical particles contributes a factor 1/j! to the annihilation rate. Note as a consequence that the units of σv n−1 depend on n; in natural units, σv n−1 has mass dimension −(3n − 4).
Suppose the annihilation has m final-state particles (here m is an index, not to be confused with the particle masses), and S f is the symmetry factor associated with identical particles in the final state, receiving a factor of j! for every set of j identical particles. Let us label the incoming momenta as p 1 · · · p n and the outgoing momenta as p n+1 · · · p n+m . Then the relation between the n → m matrix element M (as computed from Feynman diagrams) and σv n−1 is given by (e.g. [13]): (2.1) where |M| 2 is the squared matrix element averaged over the degrees of freedom in both the initial and final states, f k is the distribution function of the kth particle in the initial state, and each g i factor counts the degrees of freedom (e.g. spin) in the state labeled by i.
In the two-body case, if the initial particles are highly non-relativistic and final-state particles are highly relativistic, σv is approximately independent of the relative velocity of the annihilating particles (and hence equal to σv , which is approximately independent of the velocity distribution) for s-wave annihilation, where M is velocity-independent. Likewise, in the general case we can separate the integral into parts pertaining to the initial state and final state:

2)
Here dΠ m is the standard m-body Lorentz-invariant phase space, set by the total incoming momentum n i=1 p i ; in particular, for a 2-body final state in the center-of-momentum (COM) frame, it takes the form dΠ m |M| 2 = 1 16π 2 | p|/ √ s dΩ|M| 2 , where p is the 3-momentum of either of the two final-state particles, and s = ( n i=1 p i ) 2 . (Note that we do not include the g i factors in this phase space, but track them separately.) If the final-state particles are much less massive than the initial center-of-mass energy √ s, and hence are relativistic, and |M| 2 is momentum-independent, then dΠ m |M| 2 is approximately independent of the initial momenta when the initial particles are all non-relativistic and √ s is dominated by their masses. In this case, the integrals over the initial particle phasespace distributions can be performed trivially, yielding factors of the initial-particle number densities, and as in the 2-body case, we find that σv n−1 is independent of the characteristic velocity of the colliding particles.
If the final-state particles are not relativistic, the two-body final-state phase space will contain a factor of their 3-momentum, leading to a phase-space suppression in the annihilation rate as usual. In this case, or if |M| 2 has a non-trivial momentum dependence, then the full integral must be performed.

Thermal freezeout estimates
We now review how the standard thermal freezeout scenario is modified when the dominant interaction with the radiation bath is through a process involving n DM particles, following the estimates of Ref. [2]. For the moment, let us assume that the annihilation products are in thermal equilibrium with the SM bath. We will find it illustrative to write our results in terms of the DM mass, the temperature of matter-radiation equality (which sets the DM abundance), and the Planck mass.
Freezeout occurs when the annihilation rate is comparable to the Hubble rate. During radiation domination, the Hubble rate is given by H 2 ∼ 1 m n+1 χ , = (g * ) 1/2 f (g * ) eq (g * S ) f (g * S ) eq For the estimates contained in this section, we will assume g * S ≈ g * at all times (however, any numerical calculations track both separately); thus we obtain: σv n−1 ∼ (g * ) 3 2 −n f x −5+3n We see that in the familiar n = 2 case, the required cross-section is essentially independent of m χ , being given roughly by 1/(m Pl T eq ), and more accurately by (g * ) −1/2 f x f multiplied by this value. In the n = 3 case the required cross section scales as m −2 χ , i.e. for smaller masses a larger cross section is needed; for n = 4 the scaling is m −4 χ . If we estimate the rate as σv n−1 ∼ α n /m 3n−4 χ , on dimensional grounds (assuming a tree-level cross section and that the only relevant mass scale for the annihilation is comparable to the DM mass), we obtain an estimate for the required DM mass: For example, for 3-body annihilation, we obtain m χ /α ∼ 2 GeV (x f = 1), 40 MeV (x f = 20), 4 MeV (x f = 100), ignoring the g * factors. In general we expect perturbative 3-body annihilation to require sub-GeV DM masses if the DM is a thermal relic, and 4-body annihilation to require sub-MeV DM masses. It is also possible that the thermal annihilating component comprises only a small fraction η of the DM. In this case, the calculation proceeds as previously, except that T eq in Eq. (2.5) must be replaced with ηT eq . Thus the thermal cross section becomes: σv n−1 ∼ (g * ) 3 2 −n f x −5+3n Since the annihilation rate scales as σv n−1 (ηρ χ ) n , we see that overall the annihilation rate for a thermal relic is expected to be directly proportional to η, independent of n, in contrast to the η n scaling if the annihilation rate is held fixed and not required to be thermal.

Maximum density at late times
We now outline some generic upper limits on the signal from n-body annihilation processes. Note that the annihilation rate per particle per unit time Γ χ must be larger at freezeout than at any later time, if the mechanism for freezeout is annihilation and the DM density is not rapidly changing. In this case, at freezeout the annihilation rate is similar to H f , and since H is monotonically decreasing after inflation, H f is larger than H at any later time. Thus if the annihilation rate at any time and in any region is comparable to its value at freezeout, it will also be faster than the Hubble expansion at that time, and the DM density in that region will be rapidly depleted (in a second "freezeout"), reducing Γ χ until it is below H. If the annihilation cross section does not decrease after freezeout, it follows that the equilibrium density in any region after freezeout must always be lower than the cosmological density at freezeout. From this argument, we can set an upper limit on the (steady-state) density for the DM in any region after freezeout, ρ χ < ρ f , where ρ f is the density at freezeout, provided the annihilation cross section does not decrease at late times. In fact, from the same argument, the steady-state density at redshift z satisfies ρ χ (z) n−1 < ρ n−1 f H(z)/H(z f ). Furthermore, if the annihilation rate scales only with density, it follows that the requirement that n-body annihilation dominates over other annihilation channels at late times is always more stringent than the requirement that it dominates at freezeout (except possibly for transient periods when the DM density is being rapidly depleted). However, there can be exceptions to this latter rule if the annihilation rate does not only scale with density (e.g. if it also scales with velocity).
For fermionic DM we can set another generic constraint on the maximum DM density from phase-space considerations; this is a form of the well-known Tremaine-Gunn bound [22]. The phase-space distribution function satisfies f c ≤ 2 from the Pauli exclusion principle. Consider a region of density ρ χ , corresponding to a number density n χ = ρ χ /m χ ; then the phase-space distribution function can be approximated as f c ∼ ρ χ /(m 4 χ v 2 3/2 ) (in the nonrelativistic regime). Thus we find ρ χ m 4 χ v 2 3/2 . Note that if the two-body and n-body cross sections have "natural" scaling, i.e. σv n−1 ∼ α n /m 3n−4 χ , then it follows that the ratio of annihilation rates between the n-body and 2-body cases satisfies: (2.10) Thus we see that the natural suppression of n-body annihilation relative to 2-body, at least for fermionic DM, is a strictly stronger effect than the suppression of p-wave annihilation relative to s-wave, which corresponds to a factor of v 2 . However, this argument can break down if the n-body annihilation rate increases at low velocities, or otherwise does not satisfy the "natural" scaling estimate.
3 Low-velocity enhancements for n-body processes Above we have assumed that the annihilation is dominated by tree-level processes, and the mass scale for the effective operator is similar to the mass of the DM. However, there could be a large enhancement to the annihilation rate due to resonance effects or the presence of lighter particles in the spectrum; in the latter case, attractive interactions between the DM particles could produce a low-velocity enhancement to the annihilation rate (as discussed for Sommerfeld enhancement of two-body annihilation in e.g. [20,23,24]). Relatedly, even if two-body processes dominate during freezeout, n-body interactions might come to dominate at late times if the velocity scaling is sufficiently strong. Therefore one question we would like to address is what limits can we place on such enhanced interactions.

The unitarity bound
For two-body scattering in the non-relativistic limit, partial-wave unitarity imposes an upper bound on the cross section for annihilation in the lth partial wave, which for distinguishable particles in the initial state is given by: where k is the momentum of one of the incoming particles in the center-of-mass frame (or alternatively, the momentum of an incoming particle scattering in a stationary potential). For n > 2, a similar bound can be set for non-relativistic scattering in hyperradial potentials [25]; if k is now the magnitude of the hyperradial momentum coordinate, the 1/k 2 scaling is replaced by 1/k 3n−4 , and the prefactor is also modified. The "event rate constant" of [25], equivalent to our σv n−1 , scales with an extra power of k/µ n , where µ n is the n-body reduced mass, giving an overall momentum scaling for the lth partial wave contribution of σv n−1 l ∝ k 5−3n in the unitarity-saturating limit. The same unitarity scaling can be obtained from the optical theorem, as demonstrated in [26]. Let us focus on n → 2 processes, with an initial state denoted i and a final state denoted f . The matrix element for such a process is the same as for the inverse 2 → n process, f → i. Considering the f → f forward scattering, we can use the optical limit to set an upper bound on this matrix element: where the sum over initial states corresponds to a sum over spin configurations (and any other relevant quantum numbers), and dΠ n is the standard Lorentz-invariant n-body phase space (as in Eq. (2.3), but with the i = n + 1 . . . n + m labels replaced with i = 1 . . . n). Using this inequality and taking the non-relativistic limit to simplify Eq. (2.1) for the specific case of m = 2, one can obtain an upper bound on σv n−1 . We leave the details of the computations to Appendix A and only report the main result here: where, as before, g i counts the degrees of freedom in each state i. We see the expected scaling with T −(3n−5)/2 (∼ k −(3n−5) in the non-relativistic limit). For n = 2, in the simplified case where m i = m χ for every i, we obtain: and for n = 3, we obtain: in agreement with [26] up to the g i factors (set to 1 in that work) and the S i factor (absorbed into the definition of σv n−1 ). The 1/T 2 unitarity scaling for the rate coefficient for three-body processes has been studied in the context of ultracold bosons (e.g. [27,28] and references therein). For example, this scaling has been observed both theoretically and experimentally in dimer formation processes, where the initial state consists of three free particles and the final state contains a dimer plus a free particle [29,30]. This scaling corresponds to a very large s-wave scattering length for the two-body s-wave interactions; at sufficiently low velocity/temperature, the cross section saturates at a scale set by the scattering length.

Implications for cosmology
For a n-body DM annihilation process, let us consider the scaling of the annihilation rate with redshift before the onset of structure formation, in the regime where the unitarity limit is saturated. While the DM remains at the same temperature as the SM radiation bath, the DM temperature redshifts as (1 + z), the typical velocity redshifts as (1 + z) 1/2 , and so the unitarity scaling σv n−1 ∝ k −(3n−5) corresponds to σv n−1 ∝ (1 + z) 5 2 − 3 2 n , and the annihilation rate per particle Γ χ ∝ (1 + z) Once the DM temperature decouples from the cosmic microwave background (CMB) temperature, we instead have k ∝ 1 + z, and consequently Γ χ ∝ (1 + z) 5−3n (1 + z) 3(n−1) = (1 + z) 2 . Note that H scales as (1+z) 2 in the radiation-dominated regime and (1+z) 3/2 in the matter-dominated regime, so Γ decreases at least as rapidly as H(z) with decreasing z, and (at least once the unitarity-scaling regime is entered) there is no possibility of an early freezeout followed by a later recoupling. After kinetic decoupling of the DM, the scaling in the unitarity regime is independent of n, and thus matches that of Sommerfeld-enhanced two-body annihilation.
In general, we expect that below some saturation velocity, the cross section will fall below the unitarity limit, with a weaker velocity scaling (or no velocity scaling, or even a low-velocity suppression). If σv n−1 saturates (becoming constant) at some velocity, then below that velocity, the signal will only scale with ρ n rather than with velocity. In this case, for indirect searches probing the low-velocity regime, the effect of the velocity scaling at intermediate temperatures (above the saturation point) will be to greatly increase the normalization of the annihilation cross section σv 2 relevant for such searches, compared to the naively expected thermal value. This will be our standard assumption -that the cross section is saturated -in the indirect-detection constraints that follow. If saturation does not occur, or has not occurred for the velocities of interest for indirect detection, then the constraints we present in this work can still be used as estimates if the typical velocity of DM particles in the system of interest is known, but a detailed calculation would require inclusion of the full velocity dependence (as done for two-body annihilation in e.g. [31]).

Calculating signals from dark matter halos
A crucial ingredient in determining annihilation signatures is the contribution from low-mass halos, either isolated halos or subhalos within a larger system. This contribution is generally parametrized by the "boost factor", i.e. the ratio of the true signal to that expected from the smooth background DM density (either the cosmological DM density, for the case of probes of cosmological volumes, or the smooth density profile of the main halo). In the CDM paradigm, the smallest halos are also the oldest and most dense (having formed when the universe was denser); for few-body annihilation, the additional density scaling relative to two-body annihilation further enhances the signal from these small halos.
Unfortunately, this contribution is difficult to model precisely, as the small halos in question are below the resolution of cosmological simulations. Thus we will first describe a benchmark estimate for the boost factor, and then discuss its possible uncertainties.

Calculation of the average annihilation boost from isolated halos
In order to estimate the enhancement of the signal due to the clumpiness of the universe, we first need to obtain the expected number density of structures over a wide range of masses, for which we use the standard Press-Schechter (PS) formalism [32]. Let us start by introducing the matter power spectrum at redshift zero, which is given by: Here k is the comoving wavenumber, P 0 (k) is the primordial power spectrum, and T (k) is the standard matter transfer function taken from CAMB [33] with the fit for small scales provided by Ref. [34]. We use cosmological parameters Ω DM h 2 = 0.12, Ω b h 2 = 0.022 and H 0 = 67.27 km s −1 Mpc −1 , consistent with Planck results [35]. T χ (k) encodes the suppression of power at small scales due to the free-streaming of DM particles. T χ (k) is in general model-dependent, and in particular depends on the temperature of the DM during cosmic history. For standard WIMP-like DM particles one can compute T χ (k) analytically (see e.g. [36]). However, for the dark-sector models that we consider, we will simply approximate the suppression factor by an exponential cutoff at a characteristic wavenumber: The cutoff k c depends on the details of the kinetic decoupling between the DM and the SM.
In what follows we will study the effects of varying k c between 10h/Mpc and 10 7 h/Mpc.
Direct observations of the Lyman-α forest at k 10h/Mpc, by HIRES/MIKE [37] and XQ-100 [38], have been used to set constraints on models of warm DM (WDM) [39] and fuzzy DM (FDM) [40] via their suppression of the matter power spectrum at small scales. The WDM mass scale is currently constrained to be (conservatively) greater than 3.5 keV, whereas the conservative lower mass bound in the FDM case is 20 × 10 −22 eV (both limits are at 2σ). We note that k c ∼ k 1/2 , where k 1/2 is defined as the scale where the power drops by a factor of 2 relative to the CDM case; for the FDM model, k 1/2 is estimated to satisfy k 1/2 ≈ 4.5(m χ /10 −22 eV) 4/9 /Mpc [41]. Taking the lower mass bound quoted above, we obtain k c 17/Mpc ≈ 24h/Mpc, for the FDM model. Since we are relying on analytic estimates to translate between mass and k c , and in any case, the cutoff need not have precisely the exponential form assumed for our transfer function, we consider k c 10h/Mpc. The largest cutoff we consider, k c ∼ 10 7 h/Mpc, is comparable to the cutoff for a conventional WIMP scenario [36].
The primordial power spectrum is given by: where A s 2.21 × 10 −9 is the amplitude of the primordial power spectrum, k 0 = 0.05 Mpc −1 is the pivot scale and n s 0.969 is the primordial scalar spectral index; we have taken the Planck 2015 best-fit values for the parameters [35].
The variance of matter density perturbations for each halo mass can be obtained by smoothing the power spectrum over scales shorter than the size of the halo, i.e.: where k min ∼ 10 −4 h/Mpc corresponds to the largest observable scale, and W (kR) is a window function that describes the size of the halo. The simplest window function would be a top-hat function in Fourier space. In [42] it has been argued that this choice of window function is indeed (at least marginally) better than other possibilities when the PS formalism is compared to the simulations. Hence we only consider the sharp-k window function: where Θ is the Heaviside step function and R is a comoving length-scale associated with the halo. The disadvantage of this window function, however, is that the mass assignment to halos is not transparent, and a free parameter is needed that has to be fixed by calibration against simulation. Ref. [42] finds that the halo mass is related to the length-scale R by: with C = 2.7, where ρ m,0 is the matter energy density today. This mass assignment indicates that the range over which we change the power spectrum cutoff k c = (10 − 10 7 )hMpc −1 , roughly corresponds to halo masses in the range M c ∼ (3 × 10 −8 − 3 × 10 10 )M ; below which the structures are expected to be washed out due to the free streaming. The mass function of halos (comoving number density of halos per logarithmic mass scale at each redshift) can then be written as Here ν(M, z) = δ 2 c (z)/σ 2 (M ), where δ c (z) = 1.686/D(z) is the time-dependent critical density contrast, and D(z) is the linear growth factor normalized to 1 at the present time. For the growth factor we do the following integral numerically [43,44]: where H(z) is the Hubble parameter, H 0 is the Hubble parameter evaluated today, and Ω m is the present-day ratio of the matter energy density to the total energy density. The growth factor, in our convention, is then defined by D(z) = I(z)/I(0). For ellipsoidal collapse, the function F (ν(M, z)) is given by: with A e = 0.3222, p e = 0.3 and q e = 1 [42]. Note that for the top-hat window function in Fourier space Eq. (4.5) the halo mass function simplifies to where R(M ) is the comoving halo length-scale related to its mass via (4.6). For computing the boost factor for DM interactions, we also need some information about the halo profile. It is easy to check that the NFW profile [45] leads to a diverging integrated annihilation rate for n-body annihilation with n > 2, due to the 1/r cusp in the DM density at the center of the halo. One can modify the profile around the center by assuming a smooth core, but the actual size of the core is quite uncertain. For the majority of our studies, we will instead consider the Einasto profile [46], which describes halo profiles in N-body simulations quite well, and for which the annihilation signal is well-defined: (4.11) Here ρ −2 and r −2 are the density and radius at which the slope of the logarithmic density (d log ρ/d log r) is −2, and α e is the Einasto shape parameter. For a given halo mass, ρ −2 has to be determined by requiring that the integral of density over the halo volume gives a consistent mass. However, note that ρ −2 is a prefactor which will be cancelled out in boost factor calculations, hence it is an irrelevant quantity there. We will discuss the other Einasto parameters below.
Generalizing the results of [47], the boost factor for an individual halo in the presence of a general interaction (which need not be number-conserving), with n initial-state particles, is given by 1 : in which we have defined the following averaging of an arbitrary function G(r, M ) over the halo profile G(r, M ) = 4π with V h = 1 = 4πr 3 200 (M, z)/3. Here r 200 is the approximate boundary of the halo, the virial radius, which is defined as the radius at which the average energy density of the halo (within that radius) becomes ∆ h = 200 times larger than the total background energy density, i.e. (4.14) Note that ρ h depends on redshift through the Einasto parameters. The parameters r 200 and α e both need to be found by calibrating to either simulations or observations, the former of which is usually reported as the concentration parameter defined by c(M, z) = r 200 /r −2 . We will discuss their explicit form later but, for now, it suffices to notice that both parameters depend on the halo mass as well as redshift.
Substituting the Einasto profile (Eq. (4.11)) into Eq. (4.12), we obtain the following simple analytic expression for the halo's boost factor: 15) in which γ is the lower incomplete gamma function defined by where Γ(x) and Γ(x, y) are the gamma and the incomplete gamma functions, respectively. The universal boost factor due to the non-linear structures can then be written as: where the factor (1 + z) 3 appears because the mass function Eq. (4.7) is the comoving number density of halos. M min is determined by the cutoff in the matter transfer function k c ; using the sharp-k filter we have M min M c = 4π 3 ρ m,0 (C/k c ) 3 with C = 2.7, following Eq. (4.6). In principle we could integrate down to M min = 0 without large errors, as the exponential suppression in the mass function removes halos below M min from the integral in any case. In practice, however, we will take M min = 10 −11 M sun as a somewhat arbitrary and conservative choice. Note that the integral Eq. (4.17) quickly saturates as a function of M min so that any choice of M min (provided that M min 0.1M c ) is equally appropriate. Note that the above equations use the matter energy density, ρ m , not the DM energy density, as baryonic matter would also contribute to the gravitational potential.
The last integral in Eq. (4.17) has to be performed numerically. Note that this boost factor neglects the contribution from the smooth background in the numerator, summing only over the contributions of collapsed halos. This is valid when most annihilation occurs in halos, but at high redshifts (prior to structure formation) or in regions where most halos have been disrupted, the boost factor instead approaches unity. To capture both limits, we define the effective CDM energy density as follows: ρ eff (z) = (1 + B(z)) 1/n ρ χ (z). (4.18) Thus the n-body annihilation rate averaged over a region scales as ρ n eff , rather than ρ n . Fig. 2 shows the results for ρ eff as a function of redshift for n = 2, 3 and 4 for different parametrizations of the concentration parameter (which will be discussed below).

The effect of substructures
Within any halo there can be many subhalos that will enhance the boost factor. To estimate their impact, let us assume that all substructures are located around the virial radius of the host halo. This approximation yields a two-fold simplification: first, when integrating over the region of the halo where the subhalos are present, we can neglect the contribution from the main halo. Second, since the substructures are rather far from the center of the host halo, we can neglect tidal effects from the main halo, which are expected to truncate the substructures such that their true radius differs from naive estimates of the virialized radius [48,49].
Based on the above simplification, following the same steps that gave us the universal boost factor (Eq. (4.17)), we can write the halo's boost factor modified by substructure (B h ) as follows:B In the above relation we have assumed that ρ(r) = ∆ h ρ c (z), independent of the halo or substructure mass. M s is the subhalo mass, M sub min is the minimum halo mass set by the smallscale cutoff in the power spectrum (i.e. we set M sub min = M c = 4π 3 ρ m,0 (C/k c ) 3 as a natural choice), and M sub max is the maximum mass of subhalos (which is a function of the mass of the host halo). N s (M, M s ) is the total number of subhalos with mass larger than M s in a host halo with mass M . Ref. [50] obtains a simple fitting formula for N s (M, M s ) from simulations, which we will follow here: where µ = M s /M , and we set µ 1 = 0.01, µ c = 0.1, κ 1 = −0.94 and κ c = 1.2. Note that these parameters in principle depend on both halo mass and redshift, but these dependences -compared to other uncertainties -have only a mild effect on the boost factor. In particular, notice that since the boost factor is most sensitive to the smallest halos, changing the cutoff µ c -which only changes the number of the largest subhalos -would not significantly modify the final results. The exponential cutoff in the subhalo number also makes the above integral insensitive to the actual value of M sub max , and we simply set M sub max = 10M µ c . On the other hand, the contribution of the substructures to the boost factor is more responsive to M sub min , although the sensitivity is still mild enough that the resulting uncertainty is negligible compared to the differences amongst the range of models we will consider. For example, if we take M sub min = M c /5 rather than M sub min = M c , ρ eff is enhanced by a factor smaller than 20% in all cases (the actual value depends on the choice of n as well as the concentration parameterization).
After computing the corrected boost factorB h we insert it into the universal boost factor (4.17) simply by replacing the bare halo boost factor B h withB h .

Einasto parameters and concentration uncertainties
So far, we have used the Einasto profile and the associated free parameters (i.e. the concentration c and the shape α e ) to obtain the boost factor. We still need to discuss how the Einasto profile's parameters change as functions of halo mass and redshift. Note that for n-body annihilation the small dense halos become increasingly important as n increases, and uncertainties in the concentration of the smallest halos have a large effect on the overall signal. Estimates for the concentration parameter often require extrapolation across a wide range of halo masses, as cosmological simulations cannot resolve the smallest DM halos. In contrast, the shape parameter seems to approach to a constant value for small halos so that the uncertainties in α e are less important, although the large extrapolation is still necessary. Ref. [51] finds simple redshift-dependent fits for α e , by calibrating to simulations 2 : α e (M, z) = 0.115 + 0.0165ν (4.21) where ν(M, z) = δ 2 c (z)/σ 2 (M ) is the universality function; small halos have ν 1 and so α e is nearly independent of redshift and mass in this case. As a result of this observation and noticing that the smallest halos have a dominant contribution to the boost factor we simply fix α e to follow the parameterization given in Eq. (4.21).
In order to estimate a plausible range for annihilation signals, we now consider several alternate parameterizations for the concentration parameter that have been developed in the literature. In the following discussions, and in the upcoming plots, we sort the concentrations that we are going to study according to the resulting typical size of the boost factor, starting from the one that has the largest effect.
According to Ref. [53], observations of cluster-mass halos (with masses roughly 10 13 − 10 15 M and redshift z 0.06) suggest the following concentration-mass-redshift relation: Ref. [54] finds a somewhat different concentration based on N -body simulations, taken from redshifts z 2 and halo masses of 10 11 − 10 15 M : (4.23) Note that both these models have a power-law dependence on the halo mass, which is likely too optimistic in the sense of predicting an overly large concentration for small halos (with the degree of the overestimate being largest in the first parameterization above). Such small halos have similar collapse times over a wide range of masses, and therefore are expected to have similar natal concentrations; thus the concentration-mass relation is expected to flatten toward lower masses [55,56]. We include these power-law models for ease of comparison for earlier work, and to bracket uncertainties in the scaling of the concentration parameter for small halos.
Ref. [51] avoids this power-law extrapolation issue by expressing the concentration in terms of the universality function ν, as defined above: . For c sim1 we used Eq. (4.23) with k c = 10 7 hMpc −1 to compute the variance and universality function ν. Note that, for the smallest halos, which are the most important contributors to the boost factor, c sim2 > c sim1 at low redshifts whereas c sim1 > c sim2 at high redshifts.
Simulations discussed by Ref. [54] also show that the concentration becomes flatter as a function of mass at higher redshifts; based on this trend, Ref. [57] proposes a flat concentration relation to serve as a maximally conservative estimate: The parameterization of Ref. [58] lies between this flat relation and the universality-functionbased c sim2 estimate. We will study all the above concentration parameterizations and compare the results; Fig. 1 compares these parameterizations at z = 0 and their extrapolation at z = 50 as a function of halo mass. Note that often NFW profiles, rather than Einasto, are fitted to the observations/simulations to extract the concentration parameters; however, the NFW and Einasto profiles are quite similar except in the cores of halos. The differences are critical for annihilation signals, but not for the bulk halo shapes.
As we discussed earlier, the boost factor is dominated by the smallest halos. This fact allows us to derive an approximate scaling behaviour of the differential boost factor as a function of mass. First of all, notice that the variance in the low mass limit is almost a constant as a function of mass because a change in the halo mass, in this limit, only effectively changes the very small-scale contribution of the power spectrum in the variance integral, Eq. (4.4), which is already quite small. This, in turn, implies that ν is almost a constant in the lowmass limit. Furthermore, notice that one can parametrize the matter power spectrum at short where P is an overall amplitude, the power-law factor comes from the transfer function with τ p 2.9 and the exponential cutoff is just T 2 On the other hand, the halo boost factor Eq. (4.15), in the low mass limit, depends on halo mass only through the concentration by a power-law relation B h ∝ c 3(n−1) . This is because γ(x, y 1) ∼ Γ(x), and for small halo masses c is large and α e is small (and almost Putting all mass dependent factors together and using Eq. (4.10) for the mass function and noticing that where we have omitted the redshift dependence as well as the overall amplitude, and we have assumed that the concentration scales as a power law with respect to mass, with some index τ c , i.e. c ∝ (M/M c ) −τc . Note that all the different concentrations discussed above satisfy the condition τ c ≥ 0. Clearly, the above approximate differential boost factor peaks around M ∼ M c , showing that the boost factor is mostly sensitive to the halos with the lowest mass. In the above estimation, we neglected the substructures. As is evident in Fig. 2, the substructures would be less important in the universal boost factor. This is because, in the small mass limit, the number density of isolated halos is typically larger than that contained in larger halos, although they are of the same order of magnitude.
However, for studying individual clusters, the effect of substructures can be significant. Note that the signal from a cluster is proportional to B h when the substructures are neglected, and is proportional toB h when they are taken into account. To illustrate the significance of substructures, it is suitable to plot the ratio (B h /B h ) 1/n as a function of cluster mass which indicates how important the effect of substructures is on the effective DM energy density. See Fig. 3 where this ratio has been depicted for n = 2, 3 and 4. It is evident that the effect of substructure in the total signal is significant.

Diffuse signals
In this section we consider various indirect-detection searches for DM annihilation in halos and the intergalactic medium, and their capacity to constrain n-body annihilation. We will summarize how to estimate the size of the DM signals in each of these channels, and then discuss the resulting limits. Because we expect the signal to rise steeply at lower masses, we will focus on searches relevant to constraining sub-GeV DM, from the cosmic microwave background and from telescopes measuring X-rays and soft gamma rays.

Parametrization of the expected signal
We will characterize the strength of the indirect-detection signals by a normalization parameter: where dE/dV dt is the local rate of energy injection into SM particles (or unstable particles which decay purely to the SM) per unit volume per unit time, and ρ χ,loc is the local DM mass density. For example, for DM in the intergalactic medium far away from any structures, the local DM density could be well approximated by the cosmological average DM density, whereas for DM annihilating within a halo, the local DM density would be determined by the halo density profile, and whether the annihilation is occurring within a denser substructure. In both cases, the annihilation signal scales with the appropriate power of the local DM density, and the ratio ξ is determined by the particle physics of the theory. This parameter ξ controls the normalization of the signals to which our indirect-detection observations are sensitive, so in the following sections, we will formulate our constraints in terms of ξ.
The dependence of dE/dV dt on redshift, position etc due to variations in the DM density is explicitly canceled out in the quantity ξ. However, ξ could still possess redshift-or position-dependence; for example, if the annihilation cross section has a non-trivial velocity dependence.
Note that ξ has mass dimension 5 − 4n, so from naive dimensional analysis, we expect it to be a rapidly varying function of the DM mass for n > 2. If a O(1) fraction of the power per annihilation goes into production of SM particles, then ξ and σv n−1 are related by ξ ∼ σv n−1 m 1−n χ . The exact relationship must be calculated for each model, and depends on combinatoric factors (i.e. how many identical particles are present in the initial state) as well as the fraction of the annihilation power that goes into SM particles. In the unitarity saturating limit, as discussed in Sec. 3, the annihilation cross-section scales as σv n−1 ∝ k 5−3n . The resulting scaling of the parameter ξ is then expected to be ξ ∝ m 5−4n where v is the typical relative velocity between DM particles.

Methodology
DM annihilation during the cosmic dark ages and the epoch of reionization can modify the cosmic microwave background (CMB), and is thus tightly constrained (e.g. [59]). These constraints depend mostly on the amount of power injected in electromagnetic channels, and to a lesser degree on the spectrum of the injected particles; for decay [60] and (two-body) s-wave annihilation, they provide broadly-applicable and nearly model-independent limits on the decay lifetime / annihilation rate, for DM masses in the keV to multi-TeV range. These limits are particularly stringent for light DM, excluding the thermal relic cross section for swave annihilation for DM masses below ∼ 10 GeV (assuming there is no substantial branching fraction into neutrinos). In this section we extend these limits to 3-and 4-body annihilation. While in the case of 2-body annihilation the dominant signal comes from z ∼ 600, from Fig.  2 we can infer that signals from the epoch of structure formation will be relatively more important for n > 2.
The CMB signal is largely controlled by the rate of energy going into hydrogen ionization at redshift z from DM annihilation, which can be written in terms of the parameter ξ (Eq. (5.1)) as: Here ρ c,0 is the critical density in the present day, Ω DM is the DM mass density as a fraction of the critical density. 3 Note that the quantity (dE/dV dt) ionH is the energy injection rate per unit volume averaged over the whole universe (in contrast with (dE/dV dt) in the definition of ξ (Eq. (5.1)), which is defined locally). The function g ionH (z) captures the non-trivial rescaling of the injected energy by the boost factor (ρ eff /Ω DM ρ c ) n , and also the amount of energy deposited into ionization at redshift z, as opposed to other channels. We will also include the effects of extra excitations and heating from the secondary products of DM annihilation, which are controlled by similar rescaling functions g c (z), where c is a channel label.
The functions g c (z) are calculated, using the results of [61], from an integral over energy injection at all previous redshifts, taking into account the time needed for the injected particles to cool and lose their energy. These functions thus depend on the boost factor evaluated at all redshifts higher than the redshift of interest z. See [61] for more details; the parameter here called g c (z) is there labeled f c (z). We will follow the methodology described in Ref. [60] to determine the CMB signatures of n-body annihilation; that is, we consider injection of photons and e + e − pairs with injection energies ranging from keV to multi-TeV scales, and redshift dependence appropriate to n-body annihilation (i.e. dE/dV dt scales as ρ eff (z) n ). We compute the impact on the CMB using the CLASS public code [62], and then perform a principal component analysis (PCA) to estimate the variance in the impact on the CMB anisotropy spectrum for injections of particles at different energies. We marginalize over the standard six cosmological parameters in the PCA (see Ref. [63] for details). We then validate the limits by a Markov Chain Monte Carlo (MCMC) analysis, using the likelihoods from Planck 2015 (TT + TE + EE, low-and high-, and the lensed C ) [35], and the MontePython public code [64].
We note that limits can also be placed by requiring that the increased ionization level not overproduce the total optical depth, and that heating of the gas by DM annihilation products not violate constraints on the gas temperature [65,66]. We have checked these limits following the methodology of [66] and using CLASS to compute the modifications to the temperature and ionization histories. In particular, we require that the high-redshift optical depth to recombination from DM annihilation satisfy [67]: and that the temperature history satisfy [68][69][70]: We take the conservative approach of ignoring non-DM sources of heating, and requiring that the heating from DM does not exceed these limits.
However, under these assumptions, we find that these limits are always weaker than the CMB anisotropy bounds, so we will focus on the latter in our discussion. We note that measurements of the 21cm signal from neutral hydrogen during the cosmic dark ages could potentially set tighter bounds on the thermal history and hence on the contribution to heating from DM annihilation [71,72], but at present the only claimed detection of such a signal [73] suggests a gas temperature lower than in the standard scenario with no DM annihilation, making the extraction of limits on annihilation somewhat scenario-dependent [74].

Principal components and weighting functions
As an illustrative example to build intuition, let us first consider the signals from 3-body annihilation with two extreme concentration models, c obs and c flat , and cutoff scales of 10 hMpc −1 and 10 7 hMpc −1 . These choices span the full range of ρ eff among the structure formation models we considered. Fig. 4 shows the function g ionH (z) and its dependence on injection energy, and redshift, for these model choices. The onset of structure formation at late times spurs a steep rise in g ionH (z) at low redshifts; this effect is more pronounced when c obs is used rather than c flat , and when the cutoff momentum scale is increased from 10 hMpc −1 to 10 7 hMpc −1 , corresponding to lower-mass halos being allowed. This is expected, as these choices give rise to larger ρ eff at low redshift, as shown in Fig. 2.
The changes to the ionization history from a large energy injection (to make the effect pronounced and easily visible to the eye) are shown in Fig. 5, for the various combinations of cutoff scale and concentration model. As a benchmark for this example, we consider annihilation that produces electrons and positrons (in equal numbers) with energies of 100 MeV; in the case of a decaying non-relativistic mediator (as we will discuss for a specific model in Sec. 7), this would correspond to a mediator mass of 200 MeV.
We note that in all cases there is an appreciable cutoff-and concentration-independent modification to the ionization history at high redshift, corresponding to the signal from the smooth DM density in the epoch prior to structure formation; the choice of cutoff and concentration affects the dominant contribution from halos at z 200.
To understand the impact of changing the energy of injected particles, or injecting photons rather than electrons and positrons, we perform a principal component analysis with respect to injection energy and species. As for 2-body annihilation [75] and decay [60], we find that the first principal component consistently dominates the variance in the 3-body case, capturing over 90% of the variance; this is true for all the different concentration models and cutoffs we tested. We can thus characterize the effects of 3-body DM annihilation (for a given structure formation model) by a largely DM-model-independent pattern of perturbations to the CMB anisotropy spectrum, together with a model-dependent normalization factor. The same approach can be applied for 4-body annihilation, and likewise we find that the first principal component consistently contributes over 95% of the variance.
The first principal component (PC), evaluated for 3-body and 4-body annihilation with different concentration models but a high k c = 10 7 hMpc −1 , is shown in Fig. 6; this curve describes the significance of the changes to the CMB as a function of injection energy and species, holding ξ (and the structure formation model) fixed. The model-dependent normalization factor is given by ξ multiplied by an injection-spectrum-dependent efficiency factor, with the latter obtained by taking the integral of the spectrum against the first PC (see [60] for further discussion).
We see that the shape of the first PC differs noticeably between different concentration models; this is due to the differing redshift dependence of the signals in these cases. At low redshifts, the universe is more transparent, and the variation in deposition efficiency between injected particles of different energies is more pronounced. For example, in the 3-body case we see that the two concentration models with a weaker signal at low redshifts have very similar first PCs, as do the two models with a stronger signal at low redshifts, but the two sets of first PCs are quite different from each other. This is because in the models with smaller concentrations, the signal is dominated by high redshifts, whereas in the higherconcentration models the signal peaks near the end of the cosmic dark ages, as suggested by Fig. 5. In the 4-body case, structure formation is relatively more important, and only in the case of the c flat concentration model is the first PC very different (suggesting that only in this case does the high-redshift signal dominate). The choice of a high k c accentuates these differences; with k c = 10 hMpc −1 , the PCs for all concentration models resemble those for the c flat concentration model, as the signal is dominated by high redshifts prior to the onset of structure formation.  We see that the estimated 2σ constraint on ξ is generically strongest for injections of 100 MeV e + e − pairs (each particle has 100 MeV of kinetic energy). Cross-checking our PCA results, using the first PC only, against MontePython, for 3-body annihilation with our bracketing choices of concentration model and cutoff scale, we find fairly good agreement, with the PCA typically giving a constraint too strong by ∼ 30 − 40%. Results are given in Table I. We also test the effect of including additional PCs and find it to be small, as expected.
To indicate the importance of annihilation during different epochs, in Fig. 7 we plot the "weighting function" W(z), as derived in [61]. This function is defined so that d ln(1 + z)ξg ionH (z)W(z) governs the size of the CMB signal for arbitrary g ionH (z), in the approximation where the analysis is truncated to the first principal component. This weighting function Injection histories Table 1: 2-σ constraints on the DM annihilation parameter ξ for 3-body annihilation, for a range of cutoff scales and concentration parameter models, as estimated by PCA by only 1 PC(first column), and 3 PCs (second column) and cross-checked using MontePython (third column). Units are MeV −7 .
is extracted from the Fisher matrix governing the effect on the CMB of energy injections at different redshifts [75].
In the cases of 3-body and 4-body annihilation, where g ionH (z) generically becomes very large at low redshifts due to the effects of substructure and the small energy injection from the homogeneous dark matter component, the weighting function approximation is generally not adequate for accurate estimates of the constraints, because small errors in W(z) at low redshifts are amplified by the large g ionH (z) prefactor in the integral. However, it still provides a valid qualitative picture of which redshifts dominate the signal under different assumptions for structure formation.
As expected from previous work [60,75,76] the function peaks at redshift ∼ 600 for n = 2, and at lower redshift (z ∼ 100−200) for DM decay (n = 1). We see that for n = 4 with a cutoff of 10 7 hMpc −1 , for all but the "c flat " concentration parameter evolution, the signal is strongly peaked at low redshifts, and small halos are expected to be very important. For 3-body annihilation, on the other hand, these assumptions on the cutoff and concentration lead to an interesting bimodal structure for the weighting function, where two peaks are potentially present; one at redshift ∼ 30 driven by structure formation, and one at redshift ∼ 800 driven by the high smooth DM density at early times. From this figure, we expect the former low-redshift peak to dominate for the c obs and c sim1 prescriptions, for a cutoff scale of k c = 10 7 h/Mpc; for the c sim2 and c flat prescriptions, in contrast, the low-redshift peak becomes negligible and the high-redshift peak dominates. Smaller k c cutoffs (i.e. more depletion of small halos) enhance the high-redshift peak relative to the low-redshift peak. These results are consistent with our inferences from Fig. 6.
Repeating the constraint calculation for different cutoffs and concentration parameters, we obtain the constraints shown in (the upper panel of) Fig. 8. For the purpose of demonstrating the dependence on the structure formation parameters, we assume an injection of e + e − pairs such that each particle has 100 MeV of (kinetic) energy, corresponding to the peak of the first PC; the limits may be translated to any other energy (or to photons) by rescaling the limit by the first PC.

Constraints from the present-day universe: general points
The early universe has high DM density and the CMB provides a sensitive observational probe, but for decay and conventional 2-body annihilation, limits from dwarf galaxies (e.g. [77]), the Milky Way (e.g. [78,79]) and galaxy clusters (e.g. [80]) are often considerably more stringent. However, these bounds depend on the energy and spectrum of the annihilation products to a greater degree than the CMB limits.
As discussed in Section 2, the mass range of greatest interest for n > 2 is 1 GeV and lower, which is below the optimal energy range for the Fermi Gamma-Ray Space Telescope and well below energy thresholds for ground-based gamma-ray telescopes. For keV-scale photons, X-ray telescopes provide sensitive constraints: Chandra [81], Suzaku [82], and XMM-Newton [83] cover the range from 0.1 to 10 keV, and NuSTAR [84] is sensitive in the 3-80 keV range. INTEGRAL and EGRET data provide sensitivity to somewhat higher-energy photons, in the hard X-ray and soft gamma-ray band.
For DM below the ∼ 100 keV mass scale, the available SM annihilation channels are rather limited at low DM velocities; only neutrinos and photons are kinematically allowed. In the X-ray range, therefore, we will focus on the limits on photon lines. If the DM annihilates to relativistic intermediate particles which subsequently decay into photons, then the velocity of the intermediate particles will broaden the resulting photon spectra (to "boxes" rather than lines), likely making them somewhat less detectable. 10 MeV-GeV electrons can inverse-Compton-scatter on starlight photons (energy ∼ 1 eV), producing a spectrum of photons peaked in the keV-MeV range, and so X-ray telescopes may also have sensitivity to DM in this higher mass range.
A detailed study of X-ray sensitivity to such spectra is beyond the scope of this paper; the line limits we show should be viewed as an upper limit on the potential sensitivity (with that sensitivity being attained when the spectrum is very sharply peaked). Our principal goal here is to estimate the approximate sensitivity of various searches to multi-particle annihilation, and to compare the sensitivity of various local observations and the CMB bounds of the previous section.

Constraints from X-ray observations of the Galactic Center
Let us first consider the possible signal from the center of the Milky Way, where the DM density is expected to be large. We expect that substructure should be tidally disrupted in the Galactic Center (GC), so we do not include any boost factor from small halos. If we assume that the energy injection from the n-body annihilation is dominated by photons of energy ∼ E γ , then if the local rate of energy injection into SM particles per volume per time is given by ξρ n χ,loc , the number of photons produced per volume per time is given by ξρ n χ,loc /E γ . The integrated photon number flux from a given field of view, adjusted for the detector efficiency, is then given by: where ∆Ω = FOV dΩE accounts for the size of the field of view and the energy-independent detector efficiency E, and J is the averaged J-factor given by, where n labels the number of initial-state particles, ρ M W (r) is the DM density profile as a function of the galactocentric radius r, and is the line of sight distance from the observer, which is related to r and the angle from the GC ψ by r = R 2 + 2 − 2R cos ψ 1/2 .
The model-independent constraint on the decay rate shown in [85] fluctuates between Γ 3 × 10 −27 s −1 and 2 × 10 −29 s −1 for the 3 -79 keV energy range, with no strong trend with energy. For 10 keV photons, the limit is approximately 10 −28 s −1 , representative of this region; we convert this bound to a limiting photon flux using Eq. (5.5), yielding F GC γ ≈ 8.3 × 10 −5 cm −2 s −1 for a decay J-factor of J = 29 GeV cm −3 kpc sr −1 . Using the 3-and 4-body J-factors given above, we convert this flux limit into a limit on ξ, and display the results in Fig. 9.

Constraints from the Galactic halo in soft gamma rays
For 100 keV-GeV DM, observations of the diffuse gamma-ray emission from our Galactic halo by hard X-ray and soft gamma-ray observatories can provide stringent constraints [78].
Ref. [78] shows conservative limits on annihilation and decay by requiring that photon signals not overproduce the total observed diffuse photons by more than 2σ in any energy bin. Since no spatial information is employed, it is again straightforward to convert the limits to the 3-body and 4-body cases by simply rescaling the bounds by the J -factor within the relevant region of interest (ROI), which controls the observed photon flux. (Note that here we include only the signal from the smooth Galactic halo, and assume that the limits given in Ref. [78] are likewise dominated by the Galactic emission.) In the following observations, J -factors are quoted in units of GeV n cm −3n kpc sr −1 for n-body processes (we assume E is constant over the relevant regions of interest throughout, and so can be set to 1 in Eq. (5.6); a constant but non-unity efficiency can be absorbed into the total exposure / observation time, which cancels out in our comparison between constraints on processes with different n).

Constraints from galaxy clusters
It is also possible we can detect a signal from galaxy clusters. However, in galaxy clusters the signal is strongly dependent on the degree of small-scale substructure. Accordingly, we will show results for a range of concentration parameter estimates and cutoff scales. The integrated photon flux from n-body annihilation in a cluster, if the signal is dominated by photons of energy E γ , can be estimated as: whereB h (M ) is the boost factor given by Eq. (4.19), d is our distance from the cluster, and M is the mass of the galaxy cluster. For comparison, we also study the signal expected from a cluster with no substructure by replacingB h above with B h , given by Eq. (4.15). XMM-Newton and Chandra have observed a number of galaxy clusters and performed searches for photon lines in the few-keV range.
Here we use Perseus as an example target. As shown in [91], the flux limit observed by Suzaku is ∼ 10 −1 photons cm −2 s −1 sr −1 , for photon energies in the range from 1 keV to 10 keV, with field of view ∼ 320 arcmin 2 centered on the cluster center. Similar flux limits, corresponding to ∼ 10 −5 photons cm −2 s −1 , have been set for a 3.5 keV line in several clusters and galaxies. For example, there have been claims of a 3.5 keV line flux in the core of Perseus at a flux of 5.2 × 10 −5 photons cm −2 s −1 [92]; the Chandra ACIS-S and ACIS-I observations yielded fluxes of 1.02 × 10 −5 photons cm −2 s −1 [93]; for M31 there is a claim of a signal with a flux of 0.49 × 10 −5 photons cm −2 s −1 [93]; for Virgo there is an upper limit 0.91 × 10 −5 photons cm −2 s −1 [92]. We will not in this work seek to explain the claims of a 3.5 keV line, but they demonstrate that a line at this flux level is potentially detectable.
We use E γ =10 keV as an example to estimate the constraints on ξ. The cluster parameters for Perseus are taken to be R 200 = 1.90 Mpc, M 200 = 7.71 M , z = 0.0183 and d = 77.7Mpc (luminosity distance), following [94] (note that Ref. [94] assumes a NFW profile for Perseus, but with respect to the total mass and virial radius, we expect the difference between the Einasto and NFW profiles to be small) with ∆ h = 200 as before. The resulting bound is shown in Fig. 8 (lower panel) for a range of different choices for the substructure models. We see that in the absence of substructure, the cluster bounds are always weaker than the CMB limits, but they rapidly outpace the CMB bounds as the amount of substructure increases.

Discussion of diffuse constraints
The indirect-detection limits on 3-body and 4-body DM annihilation from the CMB, the Galactic center, the Galactic halo and galaxy clusters are summarized in Fig. 9. The late-time constraints from the CMB and galaxy clusters are highly uncertain, due to our insufficient knowledge about the non-linear structures and substructures of the universe; the major uncertainties can be quantified by the small-scale cutoff and the halo concentration, as shown in Fig. 8. These uncertainties are translated to the width of the bands in Fig. 9. For limits from the Galactic center and halo, where no substructure is assumed, the width of the bands reflects the different density profiles considered. Taking into account the effect of substructures would result in stronger constraints (i.e. pushing down the Galactic halo / center bands in Fig. 9); however, this is a complex problem, with numerical errors in simulations potentially modifying the inferred subhalo distribution in non-trivial ways (e.g. [95][96][97]). Such a study is beyond the scope of this paper.
The relation between the energy of the annihilation products and the DM mass is modeldependent in general -for example, if the annihilation produces an unstable mediator which subsequently decays, the energy of the decay products depends on the mass of the mediator, as well as the number of particles in the final state of its dominant decays -but for the purposes of this figure, we assume the annihilation products are injected with an energy equal to the DM mass. Where photon spectral information is required, we assume either that the photon energy traces the DM mass, or that (for DM masses above 1 MeV) the photon spectrum is that produced by final-state radiation from electrons and positrons; we choose whichever spectrum yields the stronger limit. Likewise, for the CMB constraints, we choose whichever channel (e + e − or photons) yields the stronger constraint, assuming that the energy of injected particles traces the DM mass. These constraints are in general not highly sensitive to the exact energy of the injected particles, so we do not expect these approximations to strongly affect the results (although models that produce X-ray photons with a similar total energy but a much broader spectrum are likely less constrained than models that yield photon lines).
In the case of Galactic center observations by NuSTAR, the constraints fluctuate rapidly with energy, as discussed above; we do not show these rapid fluctuations, but instead take the lower limit on the decay lifetime to be ∼ 10 28 s within the NuSTAR energy range, which corresponds to a fixed energy-independent limit on ξ. Including the fluctuations in the bound would modify the limits by up to one order of magnitude in either direction in very narrow energy ranges, which does not qualitatively change the results.  Figure 9: Limits on ξ for 3-body (upper panel) and 4-body (lower panel) annihilation as a function of m χ , comparing with limits from unitarity and estimates for a thermal relic scenario. The blue band shows the range of CMB constraints (estimated from principal component analysis) for different small-scale structure models, with the edges being determined by two extreme cases: c f lat concentration with k c = 10 hMpc −1 , and c obs concentration model with k c = 10 7 hMpc −1 . Below m χ = 1 MeV we use the efficiency for injection of photons, assuming the photon energy traces the DM mass; for masses above 1 MeV we use the efficiency for injection of electron-positron pairs instead, again assuming the electron/positron energy is given by the DM mass. The green band shows the estimated constraint from X-ray observations of Perseus, for the same range of small-scale structure models. The red bands show the estimated range of constraints from NuSTAR observations of the Milky Way, spanning the cored-NFW, shallow Einasto and Einasto profiles. The gray bands show the constraints from diffuse X-ray and gamma-ray observations of the Galactic halo, summarized in Sec. 5.5, for different DM profiles (Einasto, shallow-Einasto and cored-NFW). The black lines show the approximate unitarity limit from Eq. (5.9) for different typical DM velocities v. The orange line is an estimate of ξ for a thermal relic as a function of mass, taken from Eq. (2.7) setting g * = 10 and x f = 50, and ignoring combinatoric factors. Fig. 9 also shows how the limits compare to the unitarity bound for different velocities, and to the estimated size of a thermal relic signal. We see that in general, constraining thermal scenarios requires a model with more than the minimal amount of small-scale structure. However, in the presence of velocity enhancements, even the most pessimistic (and hence robust) constraints can potentially probe the parameter space for n-body annihilation. For example, the CMB signal from high redshifts can be predicted very reliably, and generically dominates the overall CMB signal in models with little small-scale structure (see Fig. 7). Future improvements in simulations of small-scale structure may help narrow down the uncertainties and allow robust constraints to be placed on scenarios that are not near the unitarity limit.
These limits in turn let us estimate the DM mass range that could potentially be probed by indirect searches. As discussed in Sec. 3, theoretical consistency requires the annihilation cross-section to satisfy the unitarity bound, Eq. (3.3) (unless there is a long-range interaction present and many partial waves contribute to the annihilation). As long as Eq. Note that for n > 2 this upper bound on mass is only weakly sensitive to the actual observational constraint ξ obs . For example, for v 2 ∼ 10 −3 (corresponding to the typical virialized velocity of particles within halos), and n = 3, an observational constraint in the range ξ obs ∼ (10 −10 − 10 10 )MeV −7 , can constrain models with m χ 1 − 1000MeV (with tighter limits having the potential to test higher masses), whereas for n = 4 and observational constraints in the range ξ obs ∼ (10 −25 − 10 5 )keV −11 we could potentially test models with m χ 3 − 10 3 keV.

Individual sources: ultracompact minihalos (UCMHs)
As we saw earlier, the n-body annihilation signal is more sensitive to small dense halos for higher n, and in general it will be more competitive with 2-body signals in regions of high density. Thus we now study the potential observability of n-body annihilation in particularly dense regions of DM that could exist in the present-day universe; we focus here on the possibility of dense ultracompact minihalos (UCMHs).

Formation and density profile of UCMHs
A UCMH is a (hypothetical) dense object that was formed at recombination or earlier; after matter-radiation equality, the UCMH develops through the accretion of DM onto the primordial seed. This contrasts with regular hierarchical structure formation where halos form at much later times. Below we briefly review the UCMH formation and the resulting density profile, following the notation and conventions of [98].
The initial seed for the UCMH formation is the DM contained within the distance scale of a mode at horizon crossing, with initial mass: where R i is the comoving radius of the mode that has just entered the horizon at the UCMH formation redshift, and H 0 and Ω DM are respectively today's Hubble parameter and the DM abundance as a fraction of the critical energy density. The mass of the resulting object remains almost unchanged until matter-radiation equality, after which structures start to grow. Since the density perturbation contrast grows linearly with scale factor we have: The accretion of the ambient DM onto the UCMH ceases to be efficient when hierarchical structure formation begins. Consequently, the final mass of the UCMH is roughly M UCMH (z ∼ 10).
The density profile of UCMHs can be analytically estimated based on a spherical collapse model for accretion of DM layers onto the primordial seed [99,100]. The self-similar solution yields the profile ρ(r) = A(z)/r 9/4 , where [98] A and ω χ = Ω DM /Ω m is the fractional DM contribution to the total matter today. This profile is consistent with simulations performed in [101,102], in contrast with earlier simulations suggesting even steeper profiles [103].
To obtain the normalization factor A(z), we need to know how the mass and radius of the UCMH evolves with redshift. Based on N-body simulations, Refs. [104,105] find the radius of the UCMH to obey: As before, since the UCMH stops growing roughly when hierarchical structure formation starts to be effective, the final radius of the UCMH is given roughly by its value at z ∼ 10. 4 Note that inserting the redshift scaling of R UCMH (z) and M UCMH (z) into Eq. (6.3) leads to a density normalization A(z) that is approximately constant with redshift.
After the UCMH has formed with this steep density profile, various physical factors (to be discussed below) are expected to reduce the density in the core region; we will model the resulting profile as a simple density cutoff: for r c < r < R UCMH 0 for R UCMH < r.

(6.5)
Here r c is the core radius, ρ cut = A(z)/r 9/4 c is the core density and R UCMH is the virialized radius of the UCMH. The core mass is M cut = 4πρ cut r 3 c /3.
Note however that more recent works [106,107] claim that the steep r −9/4 density profile should be disrupted during structure formation, or alternatively is impossible to produce from an initial Gaussian random field, with realistic minihalos being better described by shallower NFW or r −1.5 density profiles.
We will focus in this section on the optimistic case where the DM density rises steeply with decreasing r, with the self-similar r −9/4 profile, until saturating at some core radius. Where relevant, we will briefly discuss how the calculation would change for a power-law profile with a different slope. We will also test a specific example of a minihalo population with density scaling as r −1.5 at small r, using as our model the halo simulated in [107]: (r/r s ) 3/2 (1 + r/r s ) 3/2 M /kpc 3 . with r s = 10 −4 kpc. In this example, the halos are sourced by a peak in the primordial power spectrum at a scale of 7 kpc −1 . As with the r −9/4 profile, we will also truncate this profile to a flat-density core at an appropriate radius.
In addition to the uncertainties in the density profile, the UCMH abundance is not known a priori, so the results we obtain should not be considered as robust limits on the annihilation rate. Instead, we will provide estimates for the circumstances under which UCMHs could give rise to potentially observable signals, in the presence of 3-body or 4-body annihilation.

Core radius estimation
There are several distinct physical processes that may contribute to the UCMH core formation. We shall discuss the most significant ones below. In this work, we shall ignore any interplay between these effects, and estimate the actual core size as the largest cutoff radius induced by the various processes, for the r −9/4 density profile. Where relevant, we will discuss what changes to the calculation or additional ingredients are needed to compute the core sizes for other density profiles.
r v from velocity of DM The nonzero velocity of DM particles after decoupling washes out small scale perturbations as well as the inner cusp of the UCMH. To estimate the core radius due to the peculiar velocity of particles, we employ a velocity dispersion estimated from N-body simulations, as a function of comoving radius R i [104,108]: For ρ ∝ r −9/4 the velocity dispersion can then be related to the UCMH mass by: σ v (z) 0.14 1000 z + 1 From angular momentum conservation, the mean rotational velocity of particles at radius r from this velocity dispersion is: On the other hand, the Keplerian velocity of particles that orbit at the edge of a constantdensity core is given by The core radius from the peculiar velocity of particles (r v ) can then be estimated as the radius where the typical size of the peculiar velocity becomes comparable to the velocity of infalling DM particles, i.e. v rot (r v , z) v Kep (r v , z). This results in the estimate: This radius should be evaluated at the redshift of collapse, as this is when the inner density profile of the UCMH is determined; as discussed in Ref. [98], subsequent evolution primarily modifies the outer regions of the minihalo. To obtain a conservative (i.e. maximally large) estimate of the core size, we take the redshift of collapse to be at recombination, z = 1000.
For ρ ∝ r −β profile, the power 4/7 in Eq. (6.11) should be replaced by 1/ (4 − β), and the density normalization will need to be adjusted as appropriate. Eq. (6.7) depends on the initial conditions of the fluctuations prior to the collapse, and should be largely independent of the eventual halo density profile, but the relationship between the eventual UCMH mass and the comoving scale R i in Eq. (6.2) is dependent on the same formation history that affects the density profile.
For the specific case of the profile described in Eq. (6.6), from Ref. [107] we have M UCMH = 8.1 M at the simulated redshift of z = 400. Integrating the density profile of Eq. (6.6), the radial limit of integration must be R UCMH (z = 400) = 1.4 × 10 −4 kpc to obtain the correct mass. We take the collapse redshift to be z = 400 since it is the largest redshift for which we have a mass normalization; this may overestimate the core size, but we will find that even with this conservative approach, r v is very small in this scenario and other coring effects dominate. For small radii (r r s = 10 −4 kpc, the scale radius), the profile can be approximated as a r −3/2 profile with normalization: (r/r s ) 3/2 M /kpc 3 ; (6.12) this approximation is adequate for r 3 × 10 −5 kpc. By replacing A(z), R UCMH (z) and β in Eq. (6.11) we obtain r v = 3.2 × 10 −14 Mpc. This value is well inside the region where the density scales as r −3/2 , so this estimate is self-consistent.

r ann from annihilation
The DM annihilation inside the UCMH would wash out the self-similar profile, by depleting regions of sufficiently high DM density, as discussed in Sec. 2. Hence another cutoff radius is set by requiring that the core density be low enough that DM annihilations would not deplete a O(1) fraction of the DM over the age of the universe [98]: where τ UCMH is the age of the system. For UMCHs observed today, the relevant timescale is the age of the universe, τ UCMH = 13.76 Gyr. For a density profile other than ρ ∝ r −9/4 , of the form ρ ∝ r −β , the index of the bracketed term in Eq. (6.13) should be replaced by 1/β. Likewise, A(z) should be replaced by the correct normalization of the density profile. For the profile of Eq. (6.6), we use the appropriate normalization (Eq. (6.12)) for the r −3/2 profile, and check that the inferred core radius lies within the region where this approximation to the profile is valid.
Liouville's theorem tells us that the DM phase-space distribution function f c must be less than (or equal to) the peak of the initial distribution f i max . We can go beyond the earlier simple analysis for fermionic DM by specifying an initial condition at the redshift of kinetic decoupling, where the temperatures of the DM and the SM begin to diverge. Assuming the distribution is approximately Maxwellian at this epoch, we can write: Here ρ χ is the background DM energy density and T kd is the DM temperature at kinetic decoupling. f i (q = 0) = ρ χ (T kd ) m χ (2πm χ T kd ) 3/2 is then the maximum of the initial distribution. Assuming the velocity distribution of the DM in the core is non-relativistic and approximately Maxwellian, the phase space density is given by: where v 2 is the average velocity-squared of the particles in the core. Liouville's theorem then gives: If we assume that the particles in a core with radius r L are virialized, we have K = − U /2 where K and U are the kinetic energy and potential energy of the system respectively. For the kinetic energy we have K = 1 2 N cut m χ v 2 = 1 2 M cut v 2 , whereas the potential energy is given by U = − r L 0 GM (r) r dM (r) dr dr (here M (r) = 4πρ cut r 3 /3 is the mass within a sphere with radius r). Taking the integral and using the virial theorem we find: We can now use these relations to obtain a lower bound for the core radius: Provided the decoupling occurs when the DM is non-relativistic, T kd ≤ m χ , and since ρ χ (T kd ) ∝ T 3 kd , taking T kd ≈ m χ yields a lower bound on the core size of: For fermionic DM, as discussed previously, Pauli exclusion can also put a constraint on the maximum density (the so-called Tremaine-Gunn bound), and hence the core radius; the calculation is identical to that above, except that the maximum of the initial distribution ρ χ (T kd )/(m χ (2πm χ T kd ) 3/2 ) is replaced with 2. Thus the corresponding constraint can be obtained by replacing ρ χ (m χ ) by 2m χ (2πm 2 χ ) 3/2 , obtaining the limit: .

(6.20)
We find that the core radius r T G obtained by setting f i max = 2 is much smaller than r L even for a DM mass of 1 keV. For T kd smaller than m χ this ratio is even smaller. This is to be expected since the DM cosmological energy density, once the DM becomes non-relativistic, should be much smaller than m 4 χ . So we can safely ignore this constraint on the core radius in the thermally coupled scenario.
In both cases, if the density profile immediately outside the core is ρ ∝ r −β rather than r −9/4 , then the only replacement in the bounds is that the 8/15 index is replaced by 1/(3 − β/2), and the normalization of A(z) must be modified. Again, we approximate the profile of Eq. (6.6) as a r −3/2 profile with appropriate normalization, since the core lies within the inner region.

Benchmark models
We consider two benchmarks to test for sensitivity to UCMHs. On one hand, we consider a generic thermal 3-body annihilation scenario with x f = 50, T d = m χ and coupling α = 1, for both annihilation and scattering processes, corresponding to ξ = 2.5 × 10 −11 MeV −7 . In this case, the DM mass required to obtain the correct relic density is m χ ∼ 33 MeV, using the estimates of Sec. 2.
As a second benchmark, we use a parameter point for "Not Forbidden Dark Matter" (NFDM) [13], for m χ = 0.5 GeV, ε = 10 −6 , α = 4.3 and r N F DM = 1.8, corresponding to ξ = 3.4 × 10 −15 MeV −7 , which yields the correct relic abundance and is not ruled out by current experiments. We will discuss this class of models in more detail in Sec. 7. Both models, as expected for thermal relic DM with 3-body annihilations, cannot be tested by the diffuse signal searches discussed previously. For these two model choices, in Fig. 10 we compare the core sizes originating from the mechanisms discussed above. We express the size of the UCMH by the comoving radius of the corresponding mode, R i and its mass; R i is related to the UCMH mass by M UCMH by Eq. (6.4). For the profile in Eq. (6.6) the maximum core size is set by r L , with values of r L = 1.8 × 10 −9 Mpc for the generic thermal benchmark and 4.7 × 10 −11 Mpc for the NFDM benchmark.  Figure 10: UCMH core sizes induced by various physical processes for a benchmark thermal relic model with 3-body annihilation (left) and the NFDM benchmark (right) assuming density profiles ρ ∝ r −4/9 as a function of R i associated with the UCMH (defined in Eq. (6.1)) and its mass.

Point source signals from UCMHs
The integrated annihilation power in each UCMH for n-body annihilation for ρ ∝ r −4/9 is given by: , (6.21) where in the last line we have assumed that R UCMH r c . UCMHs could appear as point sources, not associated with any known sources, in observations by telescopes across a range of energies [98]. The optimal telescope will depend (as previously) on the mass range of the DM; Fermi may have sensitivity to relatively heavy DM candidates (∼ 1 GeV and above), whereas lighter DM could produce soft gamma-ray or X-ray point sources.
If the annihilation signal is dominated by photons with energies ∼ E γ , the predicted photon number flux from a single UMCH from n-body annihilation is: where d is the distance between Earth and the UCMH. By matching this flux to the sensitivity of a particular point source search, we can define d obs as the maximum distance for a UCMH to be detectable in that search. We assume that the distribution of UCMHs tracks the DM density profile of the Milky Way, which we take to be a NFW profile with parameters chosen to match [109]: c = 18, M 200 = 9.4 × 10 11 M , r s = 17.0kpc. These parameters are consistent, within uncertainties, with more recent studies of the Galactic halo, e.g. [110]. We can then determine the expected number of detectable sources, for a UCMH spectrum dominated by a particular mass scale, by integrating the total DM mass within d obs and within the field of view of the experiment, multiplying by the fraction of DM in UCMHs (which we will denote η UCMH ), and dividing the result by the UCMH mass. For experiments such as Fermi, with a large field of view where the sensitivity can vary depending on position on the sky, we can divide the sky into pixels within which the point source sensitivity (and hence d obs ) is approximately constant, repeat this analysis for each pixel, and then sum the total number of detectable sources. The predicted distribution for the actual number of observed sources from such a population will be a Poisson distribution; we will quote only the expectation values for the number of observed sources in the figures that follow, from which the probability of seeing any specific number of sources can be immediately derived.
We consider the point source sensitivity of: • The Chandra point source search [111], which has a field of view of 0. Note that for our 3-body benchmarks, the photons produced directly by annihilation will not dominantly lie in the Chandra energy band, suppressing detectability. However, these models would generically produce electrons and positrons; when O(10 − 100) MeV electrons and positrons inverse Compton scatter ∼ 1 eV starlight, the resulting spectrum should peak in the X-ray band (the photons' initial energy being increased by a factor of O((E e − /m e ) 2 ) ∼ 100 − 10 5 ). These particles can also lose energy through (for example) bremsstrahlung and synchrotron; a full analysis of these secondary photon signals would take into account all energy loss mechanisms for the electrons/positrons, and the possibility that they might propagate sufficiently far before losing all their energy that the UCMHs would no longer appear as point sources. These effects will tend to reduce the detectability of a signal from inverse Compton scattering, but their treatment is beyond the scope of this work.
Consequently, we will show default results for two benchmark models if 1 photon in this band is produced per annihilation; note that this corresponds to a rather small fraction (∼ keV/100 MeV ∼ 10 −5 ) of the DM mass energy being converted into X-ray photons. We will also demonstrate and discuss how the limit changes if the number of photons per annihilation is increased or decreased, using the NFDM benchmark as our baseline.
The Chandra point source search is more directly applicable to lower-mass benchmarks, e.g. in the 4-body case. As we will discuss below, 4-body thermal benchmark scenarios are not detectable by Chandra, even with optimistic assumptions.
• The survey mode of Fermi over a 10-year period, which views the full sky, and has a sensitivity of order a few ×10 −8 photons s −1 cm −2 for point sources producing photons at ∼ 100 MeV energies, improving to a few ×10 −9 photons s −1 cm −2 for power-law spectra extending up to higher energies.
We obtain the relevant sensitivity maps for this case using the package fermipy [112], requiring test statistic > 25 and > 10 photons per energy bin for detected sources, and using source-class data with the P8R2_V6 instrument response functions. Note that for the 3-body thermal benchmark model we consider, the DM mass is 33 MeV, and the resulting photon spectrum is expected to peak below the energy range to which Fermi has sensitivity. For the NFDM model, the annihilation we consider is mediated by a 0.9 GeV dark photon, which is between the ω and φ meson resonances. The dominant photon-producing channels are from the production and decay of mesons, and from final state radiation from electrons and muons [113,114]; the full photon spectrum is thus quite complicated to model. As an alternative to computing the full spectrum, we consider a broad spectrum with dN/dE ∝ E −2 between 30 MeV and 450 MeV. Such a spectrum yields about 10 photons per annihilation, so approximately 1 photon per annihilation if 10% of the energy goes into photons. From [114], annihilation mediated by a 1 GeV gauge boson produces on average 0.7 photons, and 5% of the energy goes into photons, similar to the result for our simplified spectrum. Therefore, by default we assume 1 photon per annihilation; we will also show the effect of changing the brightness of a source by setting different numbers of photons per annihilation.
Our estimates for the possible number of detectable sources in these two surveys are presented as a function of UCMH mass and R i in Fig. 11 for ρ ∝ r −9/4 . We plot the expected number of observable sources for the maximally visible case where η UCMH = 1; the result for other (more realistic) choices of η UCMH can be obtained by a trivial rescaling. We have also tested the halo density profile described in Eq. (6.6) for both 3-body benchmarks, and have found that in this case the expected number of observable sources is always below 1 (although such sources could potentially contribute to the Galactic diffuse background). From Fig. 11, we can see that for the 3-body case, if fractions of the DM as small as 10 −3 − 10 −2 were comprised of UCMHs with masses corresponding to R i 10 −4 Mpc scales explicitly, M UCMH 10 2 M and steep density profiles, then simple thermal relic models could give rise to potentially observable point sources in our Galaxy.
In this calculation, we consider only galactic sources because including extragalactic sources would only be relevant if massive UCMHs (about 10 10 M ) constitute a large fraction of the DM for the NFDM model we consider. This possibility is excluded by the observed distribution of DM in galaxies and dwarf galaxies. In contrast, for the 2-body thermal model more UCMHs with mass greater than 10 6 M can be observed when including extragalactic sources.
As a cross-check, we also compute the detectable number of sources from 2-body annihilation with m χ = 1 TeV and σv = 3 × 10 −26 cm 3 s −1 , assuming a sensitivity of 4 × 10 −9 photons cm −2 s −1 everywhere in the sky, which matches the assumptions of Ref. [98]. We find that in this case the sensitivity peaks for R i ∼ 10 −3 Mpc, and O(1) visible source would be predicted for η UCMH = 10 −6 , similar to the estimates of Ref. [98]. The limit presented in that work is stronger by a factor of several than one would expect from our results; the discrepancy is likely due to the bb annihilation final state assumed in that work, which produces copious photons (which is appropriate for heavy DM), whereas by default our pipeline assumes 1 photon per annihilation.
We also test a benchmark thermal relic model for 4-body annihilation; from the estimates in Sec. 2, we choose m χ = 38 keV and ξ = 4.2 × 10 −18 keV −11 . In this case we find that the number of detectable sources is below 1, even in the maximally optimistic η UCMH = 1 scenario. Examining a lower-mass (10 keV) 4-body benchmark, we find that it is likewise undetectable; the increase in ξ associated with a lower mass scale is counteracted by an increase in the core size (in this case set by r L ) and corresponding reduction in the signal. Notice that the expected signal for 3-body thermal case is below the Fermi's sensitivity. For comparison, we also show the corresponding curve for 1 TeV thermal relic DM with 2-body annihilation for Fermi. To illustrate the effect of varying brightness of sources, we set numbers of photon per annihilation for NFDM to be 10 (dashed red line), 1 (solid red line), 0.1 (dotted red line).

Early ionization and heating from UCMHs
The possible impact of annihilation within UCMHs on the early universe has been studied by Ref. [115]. We extend this analysis using the methodology of Refs. [61,116] to determine the modifications to the ionization and thermal history, and the resulting impact on the CMB. Using the experimental limits on the optical depth and IGM temperature discussed in Sec. 5.2, we find constraints that are always weaker than the CMB bounds, although again this could change with future measurements (or confirmations of current claims [117]) of the redshifted 21cm line of neutral hydrogen.
We use Eq. (6.21) to determine the redshift dependence of energy injection from UCMHs, and then perform a principal component analysis over a set of injection models corresponding to photons and e + e − pairs injected at different energies. For each of the two benchmark scenarios discussed above, we consider the resulting redshift dependence of the annihilation signal. In contrast to the analysis of Sec. 5.2, here the energy injection does not scale with the DM density raised to a high power, as the signal comes from small dense objects that are formed prior to the recombination epoch.
We have checked that the first principal component contributes over 90% of the variance shown in Fig. 12, so it is a reasonable approximation to show results for one injection energy and extrapolate to the others using the first principal component. For large R i the core radius is determined by r v , while for small R i it is determined by r L , in the two cases considered here. Both core radii are independent of redshift (in the case of r L , there is an apparent redshift dependence due to A(z), but as mentioned previously in fact A(z) is approximately constant with redshift). Since the signal mainly comes from the UCMH core, this implies that different models (2-body, 3-body and 4-body annihilation) and different choices of R i give rise to similar redshift dependences for the energy injection, and consequently the first principal component is nearly indistinguishable between these cases. Since the redshift scaling of the injected power is dominated by the evolution of the UCMH number density, proportional to (1 + z) 3 , the redshift dependence is also the same as one obtains from conventional decaying DM; the first principal component shown here is consequently very similar to that derived previously for decaying DM [60]. Similarly, the relative strength of constraints from probes of different redshifts should trace the results for decaying dark matter; in particular, confirmation of the claimed 21cm absorption signal from the EDGES experiment [117] could improve on the CMB bounds by up to two orders of magnitude [74].
In Fig. 13, we show the CMB limits on the fraction of DM that could be comprised of UCMHs for the thermal 3-body benchmark models, as well as for the 1 TeV thermal benchmark for 2-body annihilation, using the efficiency appropriate to injection of ∼ 100 MeV electrons and positrons (other annihilation channels will relax the constraints by O(1) factors). The limit for different energies and photon injections can be extrapolated by the first principal component. For NFDM and 4-body models, the CMB gives essentially no constraints for the cases we consider. We see that for these benchmark models, where there is no low-velocity enhancement, UCMHs constituting a percent-level fraction of the DM could potentially leave detectable signals in the ionization and thermal history. We note that this statement is dependent on the steep profile we have assumed for the UCMHs; we have checked that with the halo profile of Eq. (6.6), the CMB would not be able to test the scenario with η UCMH = 1 for either benchmark.  Figure 13: CMB limit on the UCMH fraction η UCMH for 2-body and 3-body annihilation for thermal models; we use the deposition efficiency for an e + e − pair with energy 100 MeV (the results for other energies/species can be obtained from Fig. 12).

Review of the model
Ref. [13] demonstrated that a simple dark photon model can have 3-body annihilations dominant during freezeout, and labeled this scenario "not forbidden dark matter" (NFDM). In this model the DM is a Dirac fermion of mass m χ , coupled to a massive dark photon A with mass m A such that 1.5m χ < m A < 2m χ . The dark Higgs responsible for breaking the dark U (1) gauge symmetry is assumed to be heavier. The dark photon kinetically mixes with the SM photon with a mixing ε 1. Thus the part of the model Lagrangian relevant to DM annihilations is [13]: where / D = / ∂ − ig / A . The annihilation ofχχ into two dark photons is Boltzmann suppressed, and the annihilation ofχχ into SM particles is ε-suppressed. Consequently, the three-body annihilation processesχχχ → A χ andχχA → A A can generically dominate DM depletion during freezeout. Ignoring the exponentially-suppressed kinematically forbidden two-body annihilation at late times (e.g.χχ → A A ), the competing two-body channel is the εsuppressedχχ → (A ) * → SM + SM − , where "SM" denotes a charged SM species.
The processχχχ → A χ can still occur at late times, as it does not require the presence of an A bath. At zero velocity, and assuming the width of the A can be neglected, the cross section forχχχ → A χ is given by: where r NFDM = m A /m χ . Let us define: 3) Then y is a slowly-varying O(1) function of r NFDM , with y(2) = 3 √ 3/4 ≈ 1.3 and y(r NFDM ) varying between 0.5 and 1.3 for r NFDM between 1.5 and 2. We can write the 3-body cross section as: where α = g 2 /4π. We see that the 3-body cross section in this model is naturally a few orders of magnitude larger than our parametric estimate of α 3 /m 5 χ , and can be tuned to be even larger by taking r NFDM close to 2. However, taking this limit also increases the DM selfinteraction rate (∝ α 2 (r 2 NFDM − 4) −2 ), which provides stringent constraints on this scenario if the χ andχ fermions constitute 100% of the DM, and the 2-body annihilation to the SM (∝ (r 2 NFDM − 4) −2 ), both of which proceed through a near-resonant s-channel A exchange in this case. The 3-body process increases more rapidly as r NFDM → 2 than either of these 2-body processes; this reflects an enhanced low-velocity scaling in the 3-body case, which we will discuss in the next section. Nonetheless, constraints from the 2-body processes can be significant, particularly if we demand a thermal history for the DM, although they may be relaxed by taking the NFDM component to constitute only a small fraction of the DM (in the case of self-interactions) or by suppressing the coupling ε.
Generally the DM mass scale preferred for this model is between 100 MeV and 1-2 GeV, due to stringent constraints on the dark photon parameter space at lighter masses, provided ε is not too small and r is not too close to 2. As a result, we generally need not concern ourselves with constraints on warm DM or extra relativistic degrees of freedom during and after Big Bang Nucleosynthesis.
Since the dominant 3-body annihilation processes in this model are χχχ → χA and the conjugate processχχχ →χA , the power going into SM final states per annihilation (assuming non-relativistic initial states) is (m 2 A + 8m 2 χ )/(6m χ ) = m χ (8 + r 2 NFDM )/6. The annihilation rate per unit volume per unit time, summing the two conjugate processes, is σv 2 (ρ tot /m χ ) 3 /8, where ρ tot = ρ χ + ρχ (we assume ρ χ = ρχ), and σv 2 describes the rate for one of the two conjugate processes, given in Eq. (7.2) (equivalently, the rate is σv n−1 × (ρ 2 χ /2)ρχ + (ρ 2 χ /2)ρ χ /m 3 χ ). Thus we obtain, for this model: In Fig. 14 we show the allowed values of ξ, and previous experimental constraints on the parameter space, for several non-fine-tuned choices of r and ε large enough to keep the dark and visible sectors kinetically coupled through freezeout; the dark-sector coupling α is chosen to yield the correct relic density. We see that in general the allowed values of ξ are 1/(10MeV) 7 or smaller, in large part because for these values of ε, a range of stringent limits forces the mass scale for the χ and A to be 10s of MeV or greater. Comparing this range for ξ to the limits of Fig. 9, we see that the NFDM parameter space for these choices of r NFDM is not constrained by indirect-detection searches, even for optimistic choices of the small-scale structure model. However, as discussed in Sec. 6, there is the potential for a visible signal in the presence of a substantial UCMH population.
It is possible that there may be more low-mass parameter space open in variations of the NFDM scenario where the mediator is not a dark photon; for example, light Dirac-fermion DM coupled to the SM via the Higgs portal may allow for open parameter space down to ∼ 1 MeV mediator masses [118], although this scenario has not been studied in the NFDM context where 3-body annihilation is included. We leave an exploration of alternate models for subsequent work.

Low-velocity enhancement in the m A → 2m χ region
In the region where r NFDM → 2, the cross sections for 2-body annihilation to SM particles, 3body annihilation and DM-DM elastic scattering can all be greatly enhanced at low velocities. Let us first neglect the width of the A and consider only the diagrams with the strongest low-velocity divergence in the limit r → 2. The leading contribution in the r → 2 limit is given by the diagrams in Fig. 15. The divergence originates from the propagator structure in the left diagram, with a similar expression for the right one obtained by the replacement k 1 ↔ k 3 . Note that in the non-relativistic limit with m A → 2m χ , the final-state particles are also non-relativistic, as the sum of masses in the initial and final states are nearly equal. We see that both the A and χ propagators go on-shell in the non-relativistic limit, as k 3 − p 1 → (m χ − m A , 0, 0, 0) = −(m χ , 0, 0, 0), and k 1 + k 2 → (2m χ , 0, 0, 0) = (m A , 0, 0, 0).
The A propagator will be regulated by the width of the A at s = (k 1 + k 2 ) 2 . Note that even if the interaction with the Standard Model is turned off completely, rendering the A stable at its pole mass, at s = (k 1 + k 2 ) 2 ≥ (2m χ ) 2 there is a contribution to the decay width from theχχ and e + e − final state. Near the pole, the corresponding imaginary part of the 1PI diagrams is given by: Note that this width goes to zero in the non-relativistic limit where the initial particle velocities go to zero; it regulates the naive 1/v 2 dependence of the propagator to 1/v. However, the same regulation does not occur for the χ propagator; the χ is absolutely stable, and while a dark-sector χ → χA vertex could regulate the propagator at sufficiently high s, in this case s = (k 1 − p 1 ) 2 ≈ m 2 χ , and this channel is closed. The creation of a nearly-on-shell stable particle induces an effective long-range interaction, which is responsible for a large enhancement at low velocities. This behavior is similar to that studied for muon colliders [119], albeit in this case the mediator is truly stable rather than just long-lived, and unlike in the muon-collider case, the singularity does not actually enter the physical region. This latter point can be proved by noting that the fermion propagator denominator is Lorentz-invariant, and if evaluated in the frame with p 1 = (m A , 0, 0, 0), becomes m A (m A −2E k 3 ); since the energy of the outgoing fermion with momentum k 3 is greater than or equal to m χ > m A /2, the denominator is always negative, and never crosses through zero for any physical choice of momenta. Likewise, for the vector propagator with denominator (k 1 + k 2 ) 2 − m 2 A , we can work in the frame where k 1 + k 2 = 0, and then (k Figure 15: Feynman diagrams providing the leading contribution to the amplitude for χχχ → χA in the low-velocity limit.
neither denominator can ever pass through zero, there is no singularity in the physical region.
Because of the absence of a physical singularity, there is no obvious necessity that the divergence be regulated, but one might still ask if the effects that regulate the singularity in the muon-collider case could soften the low-velocity behavior in this scenario. However, while beam-width effects are important in the muon-collider context [119,120], in our case the effective beam width is always very large (corresponding to astrophysical/cosmological scales) and we do not expect it to meaningfully regulate the cross section.
Thus we expect the matrix element to scale as 1/v 4 (from the two nearly-on-shell propagators) down to the point where the A width becomes significant, and then to scale as 1/v 3 until saturating due to r NFDM = 2 (there may also be a further softening of the scaling prior to saturation, if the decay width of the A becomes dominated by decays to SM particles that are not phase-space suppressed). The rate coefficient σv 2 carries an additional phasespace factor, as the final-state particles are non-relativistic, and thus we expect scaling of σv 2 ∝ v −7 ∝ T −3.5 at intermediate velocities, regulated to v −5 ∝ T −2.5 once the A width becomes relevant, and finally flattening out to a constant saturated value set by r NFDM at sufficiently low velocities.
The regulated matrix element for the χχχ → χA process, dividing by the number of degrees of freedom in both the initial and final states, and working in the non-relativistic limit with r NFDM → 2, then takes the form: where m 2 12 = (k 1 + k 2 ) 2 and P (m 12 , k 3 ) is the combination of propagators with the A propagator regularized by the decay width: and similar for P (m 23 , k 1 ). The thermally-averaged cross section can then be written as: where the phase space integrals defined in Eq. (2.3) can be simplified as: The quantities (|k * 1 |, dΩ * 1 ) describe the momentum of particle 1 in the rest frame of 1 and 2, written explicitly as and the other kinematics in the COM frame can be determined by Incidentally, while in this case the tuning between the A mass and twice the m χ mass is purely artificial, similar resonance effects could arise naturally in the context of bound states, where the bound state B has a mass slightly lighter than the sum of its constituents' masses. Replacing A with a bound state B ofχχ, the χχχ → χA process would correspond to 3-body recombination into the bound state. The strong scaling of 3-body recombination with scattering length, and the necessity for regulation to preserve unitarity, has previously been discussed in the condensed-matter literature [121]. The formation of dark matter bound states via 3-body recombination has been studied recently [122].
However, this parameter point is strongly excluded by limits on the resonantly enhanced 2-body annihilation channel χχ → e + e − from the cosmic microwave background; the predicted cross section is σv /m χ ≈ 6 × 10 −13 cm 3 /s/GeV.
More generally, the cross section for this annihilation channel is [13]: We can approximate the CMB annihilation limits by σv χχ→e + e − /m χ 4×10 −27 cm 3 /s/GeV, where we take the limit from the Planck 2015 result [35] and assume the efficiency parameter to be f eff 0.1, based on the results for annihilation to electrons in [61]. Assuming m χ m e , and that the decay of the A to e + e − is strongly suppressed by ε and hence gives a subdominant contribution to the denominator of Eq. (7.12), 5 this limit can be estimated as: On one hand suppose the second term in the denominator of Eq. (7.13) is negligible, then we have in the limit of r NFDM → 2: (7.14) On the other hand, suppose that α is large enough that the first term in the denominator is subdominant and can be ignored. The decay width for A into χχ is given by Eq. (7.6). Taking r NFDM → 2 and s 4m 2 χ (1 + p 2 /m 2 χ ) -where p is the momentum of either of the incoming particles in the COM frame -and substituting this result into the CMB limit above, we obtain: where x = m χ /T and we have set p 2 /2m χ 3T /2. Thus tuning r NFDM close to 2 requires ε to become smaller and smaller in order to evade CMB limits from 2-body annihilation.
To avoid the strong constraints from the cosmic microwave background in the presence of resonantly-enhanced annihilation to SM particles, we are led to posit a very small coupling ε between the dark and visible sectors. This also increases the lifetime of the A , but requiring that the A decays before BBN only requires a lifetime shorter than ∼ 200 seconds, which sets: Even if this condition is not satisfied, the abundance of A is generally rather small after freezeout due to its heavy mass, and thus limits from energy injection at BBN can be evaded. A much stronger condition is that ε is large enough that the dark and visible sectors remain in kinetic equilibrium throughout freezeout; this condition in general will not be satisfied if Eq. (7.14) holds, as shown in Fig. 14. If the relevant content of the dark sector is solely the A and χ,χ fields, and the dark sector is decoupled from the SM, then it will undergo a cannibalization process [123] after kinetic decoupling from the SM, as in ELDER models [26,124]. In this case the late-time relic density depends sensitively on the DM-SM scattering, and hence on ε. Alternatively, the presence of additional light degrees of freedom (which decay to SM particles prior to nucleosynthesis) in the dark sector can provide an effective dark radiation bath, which allows a standard thermal freezeout (albeit with a different number of relativistic degrees of freedom) to proceed in the dark sector. This possibility was discussed and labeled as the "secluded" scenario in Ref. [13]; it can be implemented by adding a light scalar field that couples to the DM through a higher-dimension operator. Annihilation of the DM into the light scalars is not resonantly enhanced, and consequently can easily be made unobservable at late times as well as subdominant during freezeout. The cross-section estimate in [13], requiring that this channel be unimportant during freezeout, corresponds to σv ∼ 10 −18 − 10 −12 GeV −2 ∼ 10 −35 − 10 −29 cm 3 /s for GeV-scale DM, and σv ∼ 10 −15 − 10 −8 GeV −2 ∼ 10 −32 − 10 −25 cm 3 /s for MeV-scale DM; even for MeV-scale DM, the low end of this cross-section range is unobservable in indirect detection.
In such a secluded scenario, we can obtain a second thermal relic benchmark, which we label NFDM 2 : m χ = 2 GeV, ε = 10 −11 , α = 10, r NFDM = 2 − 10 −7 gives the correct relic abundance and has ξ ≈ 6 × 10 4 MeV −7 in the saturation region, while the two-body annihilation rate evades the CMB limit. We show the evolution of ξ with DM temperature in Fig. 16. At this parameter point, the A decays with a lifetime of τ ∼ 1 seconds. This lifetime is shorter than the age of the universe at BBN, and so is unconstrained by cosmological probes, and it is also very short compared to all relevant timescales for indirect detection.
As well as the resonantly enhanced number-changing processes, χχ elastic scattering through s-channel exchange of a A also experiences resonant enhancement at low velocities, explicitly as For the NFDM 2 parameters above, σ SI m χ ≈ 43 v 100 km/s −2 cm 2 /g, to a good approximation, for 10 −7 < v/c < 0.1. This self-interaction cross section is small enough at typical cluster velocites (∼ 1000 km/s) to evade limits originating from clusters [125,126]; however, it would be very large at dwarf-galaxy scales compared to recently estimated bounds [127], and so for consistency with these limits, the NFDM component may need to be a subdominant component of the DM in this scenario (although we caution that these limits do not include baryonic effects, which could substantially modify the impact of DM self-interactions, albeit to a lesser degree in systems like dwarfs which are not baryon-dominated [128]). In a thermalorigin scenario, this should only suppress the indirect-detection signal by one power of the small fraction, as discussed in Sec. 2.
Note that in this model, it is the interplay of 2-body and 3-body annihilation that determines the relic density; in fact, at this benchmark, the 3-body annihilation rate remains fast relative to Hubble well after the asymptotic relic density is attained. The effect of the rapid 3-body annihilations is to force Y A Y χ Y 2 χ,eq /Y A ,eq ≈ Y 3 χ , where Y indicates the number density normalized to the entropy density, and the "eq" subscript indicates an equilibrium value. Since m A ≈ 2m χ , the equilibrium number density of A scales approximately as e −2mχ/T , and this requirement becomes Y 2 χ ∝ Y A , where the proportionality factor is almost constant with T . While this condition holds, the rates of χχχ → χA and its inverse process are comparable, and the net effect on both the χ and A comoving densities is small, even though both are far from their equilibrium solutions. The Boltzmann equation thus leads to nearly-constant values for the comoving number densities of both χ and A . Provided that the 3-body annihilation decouples before the A 's decay away, the late-time χ abundance is thus set by freezeout of the 2-body annihilation. This is a special case of the general argument in [13] that the freezeout of the second-slowest process controls the late-time DM density.
Finally, more extreme versions of this scenario, with a closer tuning to the resonance, could lead to extremely large self-interaction cross sections at low velocities. In such a scenario, the distribution of this strongly-interacting component within halos would likely be significantly modified from the discussion in Sec. 4; strong attractive interactions can give rise to a "gravothermal catastrophe" where the strongly-interacting component of the halo collapses to a dense cusp, potentially even forming a black hole [129]. We leave a detailed study of this scenario to later work, but note that such a mechanism could potentially enhance the expected population of steeply cusped DM minihalos discussed in Sec. 6, even if it is confirmed (as argued in [106,107]) that such UCMHs do not form and persist in standard CDM.

Assisted dark matter
A related class of models was suggested in Ref. [10], which also invokes multi-body annihilation of light DM to obtain the correct relic density, with a focus on keV-MeV scalar DM. In particular, this work claims to identify viable parameter space for keV-scale DM where 4body annihilation dominates freezeout. (They also present a model where 3-body annihilation dominates freezeout, but this scenario would be expected to have a very suppressed indirect detection signal, as the relevant 3-body annihilation channel involves a non-DM particle that is absent at late times.) This model posits a real scalar DM candidate φ and real scalar "assister" field S, with a dark-sector Z 2 × Z 2 symmetry. The SM fields are even under both Z 2 and Z 2 , whereas the φ and S fields are odd under Z 2 and Z 2 respectively. The dark-sector Lagrangian contains the operators: The Z 2 symmetry is then broken by adding an explicit interaction term between the S field and the SM, which destabilizes S; Ref. [10] uses the leptophilic operator: where l i denotes the SM leptons (and i labels the flavor). Ref. [10] focuses on the case where the coupling is flavor-dependent, with λ τ λ e , λ µ , in order to evade constraints from measurements of g − 2. This choice allows λ τ as large as 0.3 [10]. The S scalar in this case decays to two photons through a loop of τ 's, with lifetime roughly t S = 1 second for λ τ = 0.3. 6 In order to evade constraints on DM self-interactions, the authors of Ref. [10] also 6 Note that the related one-loop decay to tau neutrinos suffers a chirality suppression. choose λ φ 10 −6 . A number of 4-body annihilation channels involving φ and/or S particles in the initial state contribute to freezeout. However, for t t S , i.e. after the S scalars have decayed away, only the φφφφ → SS channel will contribute to indirect detection signals.
However, in our attempts to reproduce the results of this model, we identified an apparently significant omission in Ref. [10], to wit the neglect of the "forbidden" 2-body channel φφ → SS. For S heavier than φ this channel is kinematically forbidden at zero velocity, and consequently the rate of annihilation per volume per time experiences a Boltzmann suppression (including exponential factors only) of ∼ e −2m S /T , rather than the e −2m φ /T factor one would expect from ρ 2 φ . However, any 4-body process involving any number of φ or S particles experiences a Boltzmann suppression (in the rate per unit volume per unit time) of at least ρ 4 φ ∝ e −4m φ /T . (In both cases we are assuming the densities can be approximated by their equilibrium values up until freezeout.) In order for the φφφφ → SS process to be open, we require 2m φ > m S , and so the factor e −4m φ /T constitutes a stronger exponential suppression than e −2m S /T . Thus we generically expect the forbidden 2-body process to dominate over the open 4-body process. We have also checked numerically that including this forbidden process in the evolution equations dramatically changes the parameters required to obtain the correct relic density, indicating that it plays a dominant role in freezeout. The effect is to greatly reduce the coupling needed to obtain the correct relic density, which in turn strongly suppresses late-time signals from the 4-body annihilation.

Conclusion
We have explored the possible indirect-detection signatures of n-body DM annihilation processes that lead to the production of SM particles, either directly or through a decaying mediator. We have argued that such processes can be more strongly enhanced at low velocities than standard 2-body annihilation without violating partial-wave unitarity constraints; if the DM follows an approximately Maxwell-Boltzmann distribution, the unitarity-saturating temperature scaling is σv n−1 ∝ T (5−3n)/2 .
We have calculated the average enhancement of the annihilation rate (i.e. the "boost factor") arising from inhomogeneities in the DM density distribution for 3-and 4-body processes, as relevant to both the early universe and to large DM halos in the present day. These results may have applications beyond the case of n-body DM annihilation. For example, 2-body annihilation or scattering followed by prompt interaction of SM products with the gas could yield signals scaling as the square of the DM density multiplied by the gas density; to the degree that the gas density approximately traces the DM density, such signals could experience a similar 3-body scaling.
We have estimated general constraints on 3-and 4-body annihilation processes via their modifications to the cosmic microwave background and photon line signals from the Galactic Center, the Galactic halo, and galaxy clusters. In all cases, there are potentially large uncertainties associated with the DM density profile and the degree of small-scale structure; while these uncertainties are also present for standard 2-body annihilation, their importance is increased for n > 2 due to the enhanced scaling of the signal with density. We have provided both optimistic and conservative sensitivity estimates, based on a broad range of models for the small-scale structure; in particular, we have shown that in the most optimistic case, such constraints cannot probe unitarity-saturating cross sections for DM masses above ∼ 1 GeV in the 3-body case, or above ∼ 1 MeV in the 4-body case, unless the typical DM velocity in the system of interest is below ∼ 10 −5 c. Conversely, there are robust limits on unitarity-saturating cross sections for mass scales below ∼ 10 − 100 MeV for 3-body annihilation and below ∼ 10 − 30 keV in the case of 4-body annihilation.
The relative sensitivity of different constraints varies depending on the degree of smallscale structure assumed; unless there are very few small halos or their concentration is low, limits from the CMB and galaxy clusters tend to dominate those from the central Milky Way (where substructure is assumed to be depleted), and from the Galactic halo (where substructures are neglected). This conclusion may not hold true if the Milky Way density profile is steeper than Einasto or one takes into account substructures within the halo.
We have mapped out the redshifts that dominate the CMB signal, and demonstrated that while for 2-body annihilation the signal is very generically dominated by high redshifts (and thus insensitive to the modeling of small halos), for 3-body and 4-body annihilation the signal can readily be dominated by either low or high redshifts depending on the small-scale structure model.
An alternate (but potentially overlapping) class of scenarios occurs where the n-body annihilation determines the late-time DM density through thermal freezeout. The needed cross section scales approximately as σv n−1 ∝ m 2(2−n) χ ; thus for n > 2, there is no massindependent "thermal relic cross section" as in the 2-body case, and the sensitivity to thermal DM signals is greatly improved at low DM masses. In general, detecting thermal relic signals in the CMB or in local targets without low-velocity enhancement requires both a relatively low DM mass, below ∼ 200 keV in the 3-body case and ∼ 30 keV in the 4-body case, and an optimistic model for the degree of small-scale structure.
We have studied the "not forbidden dark matter" class of models as an example of a scenario where 3-body annihilation can dominate freezeout, and the annihilation produces particles that decay visibly in the late universe. We find that if the annihilation is not substantially enhanced at low velocities, the searches we have considered do not set any robust limits on the model parameter space that is not already excluded by other experiments; however, if there is an appreciable population (∼ 1% of the DM) of ultra-compact DM minihalos, we could potentially observe indirect signals from 3-body annihilation in the form of a new point source population and/or effects on the CMB. This class of models possesses a resonant regime where the annihilation signal is strongly enhanced at low velocities -even to a degree that exceeds the partial-wave unitarity bound, as many partial waves can contribute due to an effective long-range interaction. In this regime, we have demonstrated that it is possible to obtain visible signals in the cosmic microwave background, the Milky Way and galaxy clusters. However, in this class of scenarios such large enhancements also imply a resonant 2-body annihilation to SM particles, and a large DM-DM self-interaction cross section; the regions of parameter space that have not already been ruled out correspond to a secluded dark sector with a long-lived massive mediator to the Standard Model, where the 3-body-annihilating component may need to constitute only a small fraction of the DM in order to evade self-interaction limits. A Calculation of the unitarity limit for n → 2 annihilation In this appendix, we calculate the unitarity bound on σv n−1 starting from the optical theorem, Eq. (3.2), which we reiterate here: Note we have included the symmetry factor explicitly on the left-hand side, so the phase space needs no implicit restriction on the region of integration to avoid double-counting. Summing over the degrees of freedom in the final state as well, and then dividing by the number of degrees of freedom in both the initial and final states, we obtain: In the non-relativistic limit, and assuming a Boltzmann distribution, the phase space distribution function f i can be related to the number density (ρ i /m i ) by: Following [26], we can substitute this result into Eq. (2.1) to obtain, for n → m annihilation: where we have used the delta function to rewrite the exponential terms from the distribution functions in terms of the m final-state energies, rather than the n initial-state energies. Substituting Eq. (A.2) into this expression, specializing to the case of m = 2, and inserting the identity 1 = d 4 p 0 (2π) 4 (2π) 4 δ 4 (p 0 − p n+1 − p n+2 ), we obtain: The term on the second line is just the usual two-body phase space integral, which as previously can be rewritten as dΠ 2 = 1 16π 2 | p n+1 | COM / √ s dΩ COM , where "COM" subscripts indicate the center-of-mass frame. If we assume that the s-wave scattering dominates the elastic forward-scattering cross section (which is generic at low velocities, but need not hold if a long-range interaction is present), and the matrix element has no angular dependence, then our expression becomes: sT πT /(2 √ s)e − √ s/T , using the asymptotic expansion of the Bessel function K 1 for large argument, since we are working in the non-relativistic limit and hence can assume T √ s. Thus finally we obtain: ImM(f → f )(s).

(A.7)
In the specific case studied in [26], n = 3, m 1 = m 2 = m 3 = m 4 = m 5 = m χ , and all g i are set equal to 1 (so also there is only one relevant state f ), and consequently we obtain: Here we have approximated the non-exponential terms inside the integral by their values at s min , since in this case they converge to a non-zero value in that limit. This result is in agreement with the result derived in [26] up to the symmetry factors; we have checked with the authors of that work that they now agree with the 1/S f factor. Our cross-section is defined differently than in that reference, by a factor of S i , which accounts for the presence of the S i factor in the bound.
To set a constraint on the forward-scattering matrix element, we can consider Eq. (3.2) for the s-wave forward scattering, where the initial and final states are both set equal to f and each particle has COM 3-momentum of magnitude | k|: Consequently, we obtain: Writing the s-wave contribution to the elastic forward-scattering cross-section as σ f →f = | k| 2 , we recover the usual two-body partial-wave unitarity bound. Substituting into our equation for σv n−1 , we obtain: Again working in the non-relativistic approximation, the integral can be approximated as 2T ( n i=1 m i ) 3/2 e − n i=1 m i /T , yielding: σv n−1 ≤ 2 − 1 2 + 3 2 n T π −(3n−5)/2 S i g 4 g 5 g 1 · · · g n m 1 + · · · + m n m 1 · · · m n 3/2 . (A.12)