GeV-Scale Thermal WIMPs: Not Even Slightly Dead

Weakly Interacting Massive Particles (WIMPs) have long reigned as one of the leading classes of dark matter candidates. The observed dark matter abundance can be naturally obtained by freezeout of weak-scale dark matter annihilations in the early universe. This"thermal WIMP"scenario makes direct predictions for the total annihilation cross section that can be tested in present-day experiments. While the dark matter mass constraint can be as high as $m_\chi\gtrsim100$ GeV for particular annihilation channels, the constraint on the total cross section has not been determined. We construct the first model-independent limit on the WIMP total annihilation cross section, showing that allowed combinations of the annihilation-channel branching ratios considerably weaken the sensitivity. For thermal WIMPs with s-wave $2\rightarrow2$ annihilation to visible final states, we find the dark matter mass is only known to be $m_\chi\gtrsim20$ GeV. This is the strongest largely model-independent lower limit on the mass of thermal-relic WIMPs, together with the upper limit on the mass from the unitarity bound ($m_\chi\lesssim 100$ TeV), it defines what we call the"WIMP window". To probe the remaining mass range, we outline ways forward.


I. INTRODUCTION
A leading candidate for dark matter (DM) is a Weakly Interacting Massive Particle (WIMP) that is a thermal relic of the early universe [1,2]. For masses above ∼ 1 keV, such a particle behaves as Cold Dark Matter (CDM) [3], with dynamics governed by purely gravitational interactions. CDM is in excellent agreement with all large-scale observations of the universe, although there are some persistent discrepancies on smaller scales, where baryonic physics is also important.
The defining feature of the thermal WIMP is that its relic abundance is naturally explained by the freezeout process [4] with a weak-scale cross section, Ω χ h 2 ∼ 0.1 pb × c/ σv , where Ω χ h 2 ≈ 0.12 [5] is the DM density and σv is the thermally averaged annihilation cross section. The weak cross section would explain why nongravitational interactions of DM have so far evaded detection. In many Beyond the Standard Model (BSM) theories, there are WIMP candidates that naturally appear around the weak scale. While a simple thermal WIMP is by far not the only possibility for DM, it is a wellmotivated scenario that must be decisively tested.
Although there are strong limits on WIMP scattering from direct detection and WIMP production from colliders , these have not been shown to deliver model-independent sensitivity to generic thermal WIMP * Email: rleane@mit.edu; ORCID: 0000-0002-1287-8780 † Email: tslatyer@mit.edu; ORCID: 0000-0001-9699-9047 ‡ Email: beacom.7@osu.edu; ORCID: 0000-0002-0005-2631 § Email: chun-yu.ng@weizmann.ac.il; ORCID: 0000-0001-8016-2170 scenarios. The branching ratios, coupling types and signals are model-dependent, and so the lack of observations may just be due to such features. For example, there can be interference effects, momentum suppression, or velocity suppression, that make the direct detection and collider cross sections small even when the total annihilation cross section is not. In other words, there is no well-specified target scale for these experiments. There is for annihilation. Thus we focus on annihilation, and especially on the total cross section.
While the thermal WIMP hypothesis specifies the total cross section, there are no model-independent predictions for the annihilation branching ratios to specific final states. Thus, although there are many constraints on individual annihilation channels from various indirect detection searches, some of which even reach below the thermal WIMP cross section, a consistent and conservative interpretation of the data in the context of the generic thermal WIMP is surprisingly lacking; we aim to address this in this work.
To decisively test thermal WIMPs, the sensitivity on the total cross section considering all possible SM final states needs to be calculated. It must be tested for mass scales from ∼keV (minimum for CDM) to ∼100 TeV (unitarity) [50,51]. We ignore invisible final states as by definition, they cannot be tested with indirect detection. We also ignore neutrino final states, which are categorically more difficult to probe (though not impossible [52][53][54][55].) The most robust limits on thermal WIMP annihilation come from three sources. Below about 10 GeV, the strongest are from Planck measurements of the CMB [56], which is sensitive to the total ionizing energy injected. Electron and photon primary final states in-ject almost all their energy into electromagnetically (EM) interacting particles, while others inject at least about 25%. Due to the precision of Planck, these limits have low relative theoretical uncertainties, and are very robust. Furthermore, they do not depend on present-day annihilation rates. Planck limits apply to much lower DM masses with linear improvement down to the O(keV) scale. In addition, BBN does not allow a generic WIMP below about 10 MeV [57][58][59]. (For DM masses below the electron mass, there are also strong limits on gamma rays from annihilation in the Milky Way and beyond [60,61].) Above about 10 GeV, Fermi measurements of Dwarf Spheroidal Galaxies of the Milky Way, and the Alpha Magnetic Spectrometer (AMS) measurements of cosmicrays are the strongest robust limits. Fermi has the best sensitivity to photon-rich final states but not leptons [62,63], while AMS is most sensitive to leptons but not photons [64,65]. Particularly for AMS, there are considerable astrophysical uncertainties, but their effect can be mitigated by making conservative choices.
In this paper, we perform the first calculation of the limit on the total WIMP annihilation cross section, combining data on all kinematically allowed final states from Fermi gamma-ray observations of Dwarf Spheroidal Galaxies, AMS-02 positron flux measurements, and Planck CMB energy measurements. Together with the unitarity bound, this defines the "WIMP window," as discussed further below. In Sec. II, we discuss when the generic thermal WIMP can be considered excluded. In Sec. III, we describe our general approach to set a lower limit on the thermal WIMP mass. We then provide specific details for setting limits with Planck, Fermi and AMS in Sec. IV, V, VI respectively. We discuss our results in Sec. VII, and important progress for testing the WIMP paradigm in Sec. VIII, before concluding in Sec. IX.

II. WHEN IS A THERMAL WIMP EXCLUDED?
Searches for dark matter annihilation products have set strong limits in certain cases, requiring that the dark matter mass be m χ 100 GeV if annihilation proceeds solely to b quarks (Fermi ), τ leptons (Fermi ), or electrons (AMS). These are only upper limits, and only apply when the limit is dominated by favorable final states. How do we quantify, more generally, when the minimal thermal WIMP is excluded?
To meaningfully combine limits on all final states, we use the simple point that branching fractions of DM must add to 100%. At a particular mass, thermal-relic WIMP is excluded if, for the standard total cross section, no combination of final-state branching ratios is in accord with all constraints simultaneously. More generally, we define the limit on the total cross section as the largest value for which all constraints are satisfied.
A limit on favorable final states corresponding to m χ 100 GeV does not mean that lower masses are uninteresting. In fact, it's only for m χ 100 GeV that we then start to have the sensitivity to actually test the general thermal-WIMP hypothesis. This is when making informative statements (by combining branching ratios) starts, not ends. For higher masses, we learn nothing about a generic WIMP, only that better sensitivity is needed.
When quantifying when the WIMP hypothesis is excluded using experimental data, we must consider the following points.

A. Conservative Values for Experimental Inputs
First, in order to make a decisive statement that certain masses are excluded, we must choose conservative values for the various parameters that influence the results. While there are large uncertainties present in astrophysical searches, after conservative parameter choices, the residual uncertainties are low. This is discussed in Sec. VII B.

B. Precision of the Thermal Target Value
Second, we must ask how precisely the thermal relic cross section target is known [4]. The thermally averaged relic cross section can be expanded in partial waves, where x = m/T (x ∼ 20 at freezeout), v is the relative DM velocity, a is the velocity independent s-wave contribution, b contains the leading p-wave contribution, and c contains the leading d-wave contribution.
In the standard calculation, the uncertainty in the thermal relic cross section is very small (at the percent level) [4], arising from the uncertainty on the measurement of the matter density Ω m h 2 [5]. However, including order v 2 contributions to σv during freezeout can provide a 5 − 10% correction to the pure s-wave piece, based on the estimate that freezeout happens around x ∼ 20, i.e., v 2 ∼ 0.1. This uncertainty goes in only one direction, in that increasing the O(v 2 ) piece of σv means a smaller O(v 0 ) cross section is required to get the correct relic density. The exact relation is model-dependent.
The set of reasonable cases considered must include the vanilla thermal WIMP (s-wave 2 → 2 annihilation into visible final states) with a standard cosmological history. Thus, absent new data, the general lower limit cannot be stronger than what we calculate here. It can be weaker, which strengthens our point that GeV-scale thermal WIMPs have not yet been adequately probed.
For example, in the case of Majorana DM, the DM annihilation is helicity suppressed, scaling like (m SM /m χ ) 2 , where m SM is the mass of the DM annihilation product. For light final states, this renders the s-wave annihilation sub-dominant, and so annihilation may proceed dominantly in the p-wave. This leads to a difference in the late time annihilation cross section and the cross section at freezeout. More generally, if there are other annihilation channels that are relevant at freezeout but not for indirect detection (or vice versa), the cross section target can be affected [66]. For example, if there are multiple DM particles, there can be co-annihilation with an unstable species that has decayed by the time indirect probes become relevant. Considering the most generic scenario leads to a certain limit; considering a set of more general scenarios -which must include the generic one -can only lead to a weaker limit.

C. WIMP Signal Generation Uncertainties
Third, it is important to consider uncertainties in generating a WIMP signal. We use Pythia8.2 [67] to generate the energy spectra of stable secondary particles produced by DM annihilation, for the various SM final states. For center of mass energies ∼ 1 − 10 GeV (DM masses ∼ 0.5−5 GeV), there are hadronic resonances that may render the Pythia calculations less reliable [67,68]. While we will show Pythia-based results in this mass range, we will also present a general argument, independent of the details of the Pythia output, that Planck limits exclude the thermal relic cross section over this range (see Sec. IV A). When the final states are dominated by muons, electrons, and photons, it is possible to use analytic expressions for the energy spectra, which we use below 1 GeV.
There are also uncertainties in the modeling of stable secondary particles from the various annihilation channels (reflected in discrepancies between Pythia and other event generators [69,70]). However, the scenarios with largest uncertainty do not contribute substantially to our least-constrained combination of annihilation channels, and so we expect the effect of these uncertainties to be small. We also note that the additional radiative muon decays pointed out in Ref. [71] do not affect our results within our precision goal.

D. Choice of Statistical Significance
Fourth, the choice of statistical significance changes the upper limit on the annihilation cross section. In astrophysical searches, it is common to present limits at 95% C.L. Similarly, we present our results to 95% C.L., but show in Sec. VII B that increasing to (1 − 10 −7 ) × 100% C.L. (5σ) weakens limits by a factor of between about ∼1.5-10, depending on the DM mass. Overall, for a fixed C.L., we keep our calculations within the precision goal of 50%. That is, we neglect some uncertainties that we expect to affect the limit by less than 50%.

E. Degree of Belief
Even in light of the care we have taken to be conservative in all choices, one might not accept that a WIMP is ruled out if its maximum allowed cross section is "just below" the thermal-relic prediction. Accordingly, we later note the mass limits that result if one requires the clearance to be a factor of 2, which is arguably reasonable, or a factor of 10, which is clearly excessive. These lead to somewhat smaller lower limits on the WIMP mass, strengthening our point that the often-quoted ∼ 100 GeV for particular final states is too optimistic for a general limit.

III. GENERAL METHODOLOGY
We calculate the largely model-independent upper limit on the thermal WIMP cross section, correctly combining inputs from leading astrophysical and cosmological experiments. For Fermi and AMS, this limit is not just a linear scaling of individual channel limits from each experiment. This is because the limits are set on the gamma-ray or positron spectral energy distribution (SED) respectively. Introducing mixed final states will change the signal spectra for a given cross section, modifying the bin-by-bin energy flux, which is what determines the limit. For the CMB, while existing limits above ∼5 GeV scale linearly, in order to extend the limits to the sub-5-GeV DM mass range for general final states, we need to work in a regime where hadronic resonances are potentially important. We will present two methods for setting limits in this mass range.
All branching channels are considered where kinematically allowed (except neutrinos). This includes annihilation to electrons e, muons µ, taus τ , b-quarks b, gamma rays γ, gluons g, W -bosons W , Z-bosons Z, Higgs bosons H, and light quarks that are grouped into the channel q = u, d, s, c. We scan over all branching fractions, with 5% incremental difference for each channel, and 1 GeV incremental difference for DM mass.
For all combinations of branching ratios, gammaray and positron energy spectra are generated using Pythia8.2 [67].
To generate spectra below 10 GeV center-of-mass energies, we run Pythia with back-to-back beams with Beams:frameType=2, rather than the standard center-of-mass beam mode. We also increase the available phase space with PhaseSpace:mHatMin = 1.0. We use these spectra to derive limits using Fermi, AMS and Planck data, as described in the following sections. Once the limit on the annihilation cross section is below the thermal relic value [4], the particular DM mass is excluded. Note that while single annihilation channels, or limited combinations, have been previously fit to the multiple experiments simultaneously (e.g., [72][73][74][75][76][77][78]), a limit on the total annihilation cross section has not previously been constructed.

IV. PLANCK CMB LIMITS
Anisotropies of the CMB provide powerful insight into physical processes present during the cosmic dark ages. Any injection of ionizing particles, including those from DM annihilation, modifies the ionization history of hydrogen and helium gas, perturbing CMB anisotropies. Measurements of these anisotropies therefore provide robust constraints on production of ionizing particles from DM annihilation products. The most sensitive measurements to date are by Planck [56], superseding earlier measurements by WMAP [79].

A. Energy Injection from Annihilating DM
The power deposited by DM annihilation, controlled by the parameter determines the strength of the CMB limit. Here σv is the thermally averaged DM annihilation cross section and m χ is the DM mass. We calculate the weighted efficiency factor f eff by integrating our electron/positron and photon energy spectra from Pythia over the f e ± ,γ eff (E) curves calculated in Ref. [80], (3) Following Ref. [80], we neglect the contribution to energy deposition from protons and antiprotons; generally only a small fraction of the total energy of the annihilation products goes into pp production, and protons and antiprotons also deposit energy less efficiently than electrons, positrons, and photons [81]. Including these contributions would slightly strengthen the constraints. From Planck data, the 95% C.L. limit on p ann is [56] Figure 1 shows the single-channel limits on the cross section from the CMB. Below 5 GeV DM mass, as there is extra uncertainty in the Pythia spectra, we also present arguments for the thermal WIMP exclusion based on generic arguments about the efficiency and energy injection rate, as discussed below. Figure 2 shows the fraction of power proceeding into EM channels (electrons, positrons, and photons) is quite stable as a function of DM mass, and is 26% or higher for all DM masses between 5 GeV and 10 TeV and all (nonneutrino) final states. (Note these are the final branching ratios after all the unstable particles decay, not the direct branching ratios into electrons, positrons, and photons.)

B. Energy Injection Fractions
For annihilation to electrons and muons, this statement is fairly trivial, as the fraction is 100% for electrons or simply determined by 3-body kinematics for muons (resulting in roughly 1/3 of the energy going into electrons and the other 2/3 into neutrinos); the only subtlety is at high masses, where electroweak corrections allow neutrinos to be produced even for the e + e − final state. For other final states with hadronic decays, the decays will typically proceed through either neutral or charged pions. Neutral pions decay to photons with a nearly 100% branching ratio, while charged pions decay first to a muon and antineutrino (or antimuon and neutrino), with the muon then decaying to an electron + neutrino + antineutrino. In the latter process, the muon receives ∼ 80% of the pion rest energy, and then the electron carries away roughly 1/3 of the muon energy, so ∼ 25% of the pion's energy is carried by the final-state electron.
From these simple arguments, we would expect the fraction of power into EM channels to vary between ∼ 25% and 100% for non-neutrino final states, in agreement with Fig. 2. There is no reason to expect this argument to break down for DM masses below 5 GeV, although the branching ratio into hadronic versus leptonic final states may change rapidly when the DM mass becomes close to a hadronic resonance. For hadronic final states, we furthermore expect that the energy of the produced photons/electrons will peak no lower than a O(1) fraction of the pion mass (and the QCD scale); likewise, muon decays will typically produce electrons with O(10-100) MeV energies.
We therefore robustly expect that for DM masses between ∼ 100 MeV and 5 GeV, at least 25% of the DM rest energy should go into producing photons, electrons and positrons with energies above 5 MeV. Comparing with the dashed lines in Fig. 2, which are results from Pythia, confirms this argument in the hadronic resonance region. However as Pythia carries extra uncertainty in this regime, we use this estimate as a conservative cross check to set a robust constraint on light DM annihilation.
For electron/positron/photon energies above 5 MeV, the minimum value of f e ± ,γ eff is 0.32. Thus we expect f eff for any 2-body SM final state other than neutrinos to exceed f min = 0.25 × 0.32 ≈ 0.08 for DM masses in the 100 MeV -5 GeV window (for masses below this window, the results for direct annihilation to electrons/positrons/photons can be used). (As a cross-check that this is conservative, the minimum f eff value for DM masses above 5 GeV is 0.12 for the same set of channels; realistically all the electrons/positrons/photons will not be concentrated at the energies that minimize f eff .) Taking the Planck 95% confidence limit in Eq. (4), this conservative f min model-independently implies for DM masses below 5 GeV, which definitively excludes the s-wave thermal relic cross section in this mass range.
V. FERMI-LAT DWARF SPHEROIDAL GAMMA-RAY LIMITS Dwarf Spheroidal Galaxies (dSphs) of the Milky Way are one of the best DM signal targets, as according to kinematic data they are DM dense with low background. Fermi has set limits on gamma-ray fluxes from discovered dSphs [62,63], with no conclusive DM signal. DSphs provide some of the strongest limits on DM annihilating to any photon-rich final states, such as gamma-ray lines or hadronic final states.

A. Fit to data
To set limits on photons from mixed final states, we follow the official Fermi analysis on Pass 8 LAT data [63] and consider a total of 41 dwarf galaxies, both kinematically confirmed and likely galaxies 1 . Where provided, we use the measured J-factor and uncertainty. This is for 19 dwarf galaxies: Bootes I, Canes Venatici I, Canes Venatici II, Carina, Coma Berenices, Draco, Fornax, Hercules, Leo I, Leo II, Leo IV, Leo V, Reticulum II, Sculptor, Segue 1, Sextans, Ursa Major I, Ursa Major II, and Ursa Minor. For the remaining 22 galaxies not named above, spectroscopic J-factors are unavailable, and following Ref. [63] we use the predicted Jfactors with a nominal uncertainty of 0. 6  Note that four of these galaxies (Reticulum II, Tucana III, Tucana IV, Indus II) have shown a local ∼ 2σ excess in gamma rays [63], which may be attributed to DM. However, this is not globally significant, so we do not fit these excesses to a DM signal, and instead treat the measurements as exclusion bounds.
For each of these dwarf galaxies, Fermi provides the likelihood curves as a function of the integrated energy flux, where J i is the J-factor for dwarf i. Following Ref. [62], we treat the energy bins as independent, and obtain the full likelihood L i (µ|D i ), which is a function of the model parameters µ and data D i , by multiplying the likelihoods for each for the 41 dwarfs together. The uncertainty in the J-factor is included as a nuisance parameter on the global likelihood, modifying the likelihood, × 1 ln(10)J i √ 2πσ i e −(log 10 (Ji)−log 10 (Ji)) 2 /2σ 2 i as per the profile likelihood method [82]. For log 10 (J i ) and σ i , we use the values provided in [63] for a Navarro-Frenk-White profile [83]. The likelihood is maximized to produce an upper limit on the annihilation cross section at 95% C.L. Figure 3 shows our limits for the 100% branching fraction scenario. For the channels for which results are shown in Ref. [63], our results are comparable. Fermi gamma-ray line searches prohibit a large branching into gamma rays [84,85]. Note that for TeV DM masses, limits from Imaging Atmospheric Cherenkov Telescopes (IACT) such as the High Energy Stereoscopic System (H.E.S.S.) are stronger [86][87][88][89], but do not probe the thermal relic cross section. H.E.S.S. gamma-ray line searches take over where Fermi loses gamma-ray line sensitivity, and prohibits a large branching into gamma rays. Thermal relic cross section is the black dashed line [4].

VI. AMS-02 POSITRON FLUX LIMITS
AMS measurements of electron and positron fluxes and fractions [64,65] provide the strongest constraints for electron and muon final states. We use the positron flux data to set a limit on DM annihilation to all final states. As we aim to set a robust exclusion on the WIMP annihilation cross section, and cosmic-ray propagation is not precisely understood, we must take conservative parameter values at every step.

A. Cosmic-Ray Propagation
To propagate our positron spectra generated with Pythia, we use the cosmic-ray propagation program Dragon [90,91].
The evolution of the number density N i of injected electrons and positrons is given by the diffusion equation, where the convection and diffusion in momentum space have been set to zero, as it does not largely affect the spectrum in the energy region of interest [92]. In Eq. (8), D is the spatial diffusion coefficient, parametrized as where ρ = p/(Ze) is the rigidity of the charged particle with Z = 1 for electrons and positrons. The diffusion is normalized by D 0 at the rigidity ρ 0 = 4 GV. We assume the diffusion zone is axisymmetric with thickness 2z t . In Eq. (8),ṗ accounts for the energy loss; Q i is the source term, where the DM source contribution is The effect of nuclei scattering with the gas is modeled by the second line in Eq. (8).
The injected positrons are propagated using the model of Ref. [93]. This sets z t = 4 kpc, D 0 = 2.7 × 10 28 cm 2 /s, δ = 0.6, but we take the local DM density to be the conservative ρ = 0.25 GeV/cm 3 , with an NFW profile. We set the magnetic field at the Sun to be B = 8.9 µG, which means that the local radiation field and magnetic field energy density is 3.1 eV/cm 3 . As this is even higher than the conservative value of 2.6 eV/cm 3 [94], it leads to a higher energy loss rate for the positrons, which is the second-leading effect for CR propagation, behind the leading effect of the local DM density. As such, different choices of the other propagation parameters do not appreciably change the results.
The most substantial energy-loss for charged cosmic rays below about 10 GeV is due to solar modulation. The largest measured value of the solar modulation potential during AMS's data taking period of Φ = 0.6 GV is taken [95], and we employ the force-field approximation, which is valid for positron fluxes [95,96].

B. Fit to Data
We assume the data measured by AMS do not have any DM source contributions. I.e., we do not assume the additional smooth ingredient over the astrophysical background measured by AMS is from DM, as the source of the additional ingredient is unknown. As such, rather than model backgrounds, we parameterize the total AMS measurements with a degree 6 polynomial function of variable log(energy) fit to log(flux). To set the limit we perform a likelihood ratio test, where the likelihood function is where θ = {θ 1 , θ 2 , ..., θ n } are parameters in the best fit polynomial function, and the χ 2 (θ) is given by where f th i is the prediction from the modeled background (the polynomial function), f data i is the central flux energy bin of the AMS data, and σ i is the uncertainty on the particular flux value. We sum over all the AMS energy bins measured, and add the systematic and statistical uncertainty in quadrature. We then allow the parameters of the function to float within 30% of their best fit values without DM (increasing the allowed values does not weaken constraints), that allows the function to absorb a DM signal if it is preferred over the additional smooth component. We increase the DM signal normalization until the functional fit of the background plus signal to the data produces a χ 2 which has increased by 2.71, i.e., This produces an 95% C.L. upper limit on the DM annihilation channel. Figure 4 shows our limits for the case of 100% branching fractions. We check that we can reproduce comparable results to similarly conservative scenarios from Ref. [96,97]. Our results are more conservative than the weakest region of the bounds presented in Ref. [94], which arises from taking the choices for propagation parameters that lead to the weakest limit (see Sec. VII B for more details).

VII. RESULTS AND DISCUSSION
We now present our results combining all limits from Planck (Sec. IV), Fermi (Sec. V) and AMS (Sec. VI), using the method of combining all kinematically allowed branching fractions as described in Sec. III. Figure 5 shows our calculation of the largely modelindependent upper limit on the WIMP total annihilation cross section. We find the model-independent lower limit on the DM mass in the generic scenario is 20 GeV -this is where the total cross section limit crosses the thermal relic line. Following the discussion in Sec. II E, note that enforcing a clearance factor of 2 leads to a lower limit on the DM mass of about 6 GeV, and a clearance factor of 10 gives a lower limit of about 2 GeV. Also shown for comparison are the standard 100% scenarios commonly considered in the literature, τ final states probed by Fermi, and e final states probed by AMS. It is clear that the general approach of comparing favorable single-channel limits with the thermal relic cross section badly overstates the degree to which thermal WIMPs have been probed.

A. Limit on the Total Annihilation Cross Section
At different DM masses, different experiments dominate the total limit. These regions are as follows: m χ 135 MeV: In this region, the CMB is most constraining. The only available final states are muons, electrons, and photons. Therefore, we use the analytic spectra to obtain a limit, taking the least constraining of the three final states to find the combined limit. This region is excluded for generic WIMPs. The drop in the limit at ∼ 105 MeV, the muon mass, is the kinematic cutoff from a limit on muon final states, to electron and photon final states.
0.135 m χ 5.1 GeV: In this region, the CMB is most constraining. For these masses, hadronic resonances introduce extra uncertainty in spectra generated with Pythia. However, taking the conservative limit from Eq. (5), the bound still remains below the relic line. The limit shown in this region is generated using Pythia. This region is excluded for generic WIMPs.

m χ 7 GeV:
In this region, the CMB is most constraining. There are no hadronic resonances, and this limit simply comes from taking the least constrained final state from the CMB, taus, shown in Fig. 1. As shown in Fig. 5, this region is below the thermal relic line, and so is excluded for generic WIMPs. m χ 1000 GeV: In this region, the combination limit from Fermi and AMS is most constraining. This is shown as "Fermi+AMS" in Fig. 5 (note the H.E.S.S. gamma-ray line search in this region prohibits large branching into gamma-ray line photons). This region is where the total cross section limit crosses the thermal relic line, giving a model-independent lower limit on the DM mass of m χ 20 GeV.  Figure 6 shows the threshold branching fractions for fixed masses, with DM masses binned with width 1 GeV. For the lower DM mass limit of 20 GeV, these are 60% to muons, 30% to gluons, 10% to taus, and 0% to the remaining final states. Note that the exclusion branching fraction combination is not necessarily unique. For many masses there are permutations of more than one final state that give a comparable limit (i.e., swap the gluon region with b-quarks, or some linear combination of bquarks and gluons). The muon contribution however is generally present in all combinations, it is the smaller remaining combinations that vary. Muons are the least constrained among all the visible annihilation products.

B. Quantifying "Conservative"
Our limits are conservative; that is, we have consistently made choices that lead to a weaker final limit, within the parameter uncertainties. This is required to claim a meaningful robust lower limit on the DM mass. In this subsection, we detail steps taken to ensure a conservative limit for each experiment, and discuss variations from astrophysical uncertainties and modeling. Figure 7 compares our AMS limit for the χχ → ee annihilation channel with the limit obtained in Ref. [94]. This shows our limit is more conservative -i.e. weakercompared to existing bounds in the literature. The differ- ence is because we have chosen parameter values (within the allowed ranges) that lead to weaker constraints, to ensure that any resulting exclusion will be robust. To reiterate and summarize our earlier discussion, in obtaining our limit we have made the following choices:

AMS
• We do not model the backgrounds, and instead parameterize the total AMS measurements with a polynomial function. This does not require any assumptions about the interplay of background and signal propagation and their origins; we can take only the signal propagation as conservative. We check that this method produces a limit equally as weakly constraining as other works that do model both backgrounds and signals in a conservative manner [94].
• We employ high energy losses by taking a conservative choice for magnetic fields, of B = 8.9 µG at the Sun.
• We take the largest value of the solar modulation potential, Φ = 0.6 GV, measured for AMS during its data-taking period [95]. In Fig. 7, the difference in our limit at low masses and other references is sourced by this choice.
• The local DM density is constrained to the range ρ = [ 0.25, 0.7 ] GeV/cm 3 [98]. We take the low- est density of ρ = 0.25 GeV/cm 3 . This has the most dramatic impact on the limit, and is shown in Fig. 7. Other choices such as propagation model, or choice of DM halo profile, have a sub-dominant effect on our result.
These choices lead to a robust and conservative constraint within the framework of the standard diffusive propagation scenario for cosmic rays. If the propagation of electrons and positrons is substantially different than expected (e.g., [99,100]), this could impact our constraints; however, any such modification would need to obey stringent constraints from the wide range of cosmicray measurements at Earth (e.g., [101,102]). Furthermore, our modeling of the unknown background with a smooth function should accommodate many modifications to the secondary flux of cosmic-ray positrons from astrophysical sources.
AMS also reports strong limits on DM annihilating to b-quarks from their antiproton dataset. These limits are stronger than Fermi at low DM masses ( 50 GeV) [103]. We do not include the antiproton data in our combined limit. This is because the b-quark channel is not one of the key threshold channels; the weakest channels from each experiment are what sets the combined limit. Fig. 6 shows that the threshold does not contain b-quarks, and so introducing an additional analysis with stronger b-quark limits would not affect the limit on the total annihilation cross section.
Lastly, note that where AMS begins to lose sensitivity at low DM masses due to solar modulation, Voyager 1 begins to provide stronger limits [104][105][106], as Voyager 1 crossed the heliopause during data taking. We do not include Voyager 1 limits in our analysis, as the CMB is more constraining in this region.
As the AMS dataset has the largest uncertainties of all those considered in this paper, the choices discussed above have the greatest impact on the total annihilation cross section limit. Uncertainties from Fermi and the CMB are comparably negligible. Regardless, we now discuss their relevant sources of uncertainty.

Fermi
The largest uncertainty from Fermi is from the values of the J-factors. Taking a larger J-factor uncertainty of 0.8 dex does not change the Fermi limit by more than about 10% [62]. The choice of DM halo profile leads to a negligible change, as the innermost region of the halo, where the density is most uncertain, does not dominate the limits.
Note that one of the dwarfs in the nominal sample, Tucana III, shows evidence of tidal stripping. This is a likely explanation of the reported excess of gamma rays over background, rather than DM. Excluding these systems strengthens the limit, but not substantially. Using the 2015 Fermi analysis [62] instead, which did not include any systems with excesses and only kinematically confirmed galaxies, increases the lower limit on the WIMP mass by only ∼few GeV. As this is clearly not a large effect, and to be most conservative, we choose to use the most recent full dataset which gives a weaker limit due to these excesses. (Note this low significance excess can be meaningfully combined with the AMS antiproton excess [74][75][76][77][78].) The case of DM annihilation in the Milky Way halo was considered recently in Ref. [107], where stronger bounds were found for DM annihilating to b-quarks, compared to dSphs. However, similar to the AMS antiproton bounds, as Fig. 6 shows that the threshold branching combination does not contain b-quarks, introducing an additional analysis with stronger b-quark limits would not affect the limit on the total annihilation cross section.

Planck
The CMB bounds are the most robust and come with little theoretical uncertainty, especially as they do not depend on late time annihilation rates. In generating our CMB limits, the largest potential source of uncertainty comes from the energy spectra generated with Pythia. However, in the less certain hadronic resonance regime we present arguments based on ionizing energy injection fractions in Fig. 2, and using Eq. (5) confirm the limit is still below the thermal-relic line in this regime. Figure 8 shows the variation of the AMS and Fermi muon limits with higher statistical significance. We show the C.L. we take for our combined limit, 95% C.L. (2σ), as well as 99.7 % C.L. (3σ) and (1-10 −7 ) × 100% C.L. (5σ) values.

C. Implications for WIMP models
We have studied a generic WIMP: an s-wave thermal relic with 2 → 2 annihilation to visible final states, with a standard thermal history and radiation dominated early universe. While this is a largely model-independent approach, the results have important implications for DM models.

Model Types
Our approach covers models that have suppressed collider or direct detection signals. Important examples include scattering rates which are suppressed by powers of velocity or momentum, or cancellation between scattering diagrams [108][109][110][111]. For Wino or Higgsino DM candidates, scattering occurs through suppressed loops [112], but its annihilation is not suppressed.
The threshold branching fractions in Fig. 6 show that muons are the least constrained final state among all the visible annihilation products, and Fig. 5 shows the combined limit is closest to following the AMS muon limit line for masses above 10 GeV, and the Fermi muon limit line for masses above about 100 GeV. Large couplings to muons are possible within leptophilic DM models [113][114][115][116][117][118][119]. Interestingly, even if DM does not couple to hadrons at tree-level, such interactions are induced at loop level, leading to hadronic contamination of the energy spectra [119]. We do not include such effects in our spectra, but as these considerations are fairly generic, they could be used to make general statements about the WIMP.

2 → n Processes
Compared to 2 → 2 processes, 2 → n processes have energy spectra that are softened [97,120]. Therefore, such complications lead to even weaker lower limits on the DM mass compared the 2 → 2 case. This further supports the arguments in this work.
Our framework can be extended to hidden sector models with small DM-SM couplings, leading to an on-shell mediator 2 → 4n scenario, χχ → (2n × Y ) → (4n × SM), where n is the number of cascade decays and Y is a dark sector mediator. However, such a limit is only generically made assuming either the same mass mediator at each dark cascade step, or that the masses are much : Bounds on the generic thermal WIMP window (s-wave 2 → 2 annihilation, standard cosmological history), assuming WIMP DM is 100% of the DM. Shown is the conservative bound calculated in this work from data (Visibles), and the unitarity bound [50]. The remaining WIMP window is the orange line, and the white space is unprobed. Thermal relic cross section is the dashed line [4].
lower than the mass of their progenitor particle; otherwise the portion of DM energy split into each mediator's final states will be unequal [121,122], introducing extra model dependence to the calculation.

Invisibles and Sub-Dominant Density
When the limit on the total cross section is below the thermal-relic prediction, the WIMP is nominally excluded. There are two other possible interpretations.
First, the fraction below the limit can be interpreted as the fraction required to proceed to invisible final states.
Second, the strength of the limit below the relic line can also be used to set a bound on sub-dominant WIMP content. For standard indirect detection analyses for WIMP DM, the annihilation cross section and the density are often considered as independent, and are related to the astrophysical flux F as Shown is the conservative bound calculated in this work from data (Visibles), and the unitarity bound [50]. Thermal relic cross section is the dashed line [4].
where ρ χ is the DM density, and is the line of sight. The upper limit is obtained from upper limits on F , i.e., For sub-dominant WIMP DM, if the WIMP density is completely determined by the annihilation cross section, they are no longer independent, as where σv χ ∼ 3 × 10 −26 cm 3 /s is the thermal relic cross section. The annihilation flux from the sub-dominant WIMP is then Therefore, an upper limit on the flux implies σv χ 2 σv WIMP < σv limit , which provides a lower limit on the sub-dominant WIMP cross section, Figure 9 illustrates the allowed mass range for a thermal WIMP, between our new general bound from Visibles (all SM states but neutrinos) and the bound from unitarity [50], assuming WIMPs are 100% of the DM. This window is 20 m χ 1000 GeV. Cross sections below the thermal-relic line are shaded as "Overabundance", as WIMPs with smaller cross sections produce more DM than observed, which is constrained with high accuracy [5]. The unitarity upper bound at larger masses can be escaped in non-standard scenarios, such as composite DM or in the presence of extra degrees of freedom, see i.e., Refs. [139][140][141]. (Also note that in the presence of light mediators, contributions of higher partial waves to the cross section can weaken the unitarity constraint [142].) Note the bound shown in Fig. 5 has been extended here -for DM masses m χ TeV, CMB bounds are strongest. Figure 10 illustrates the allowed mass range for a subdominant WIMP, with a lower bound on the cross-section of WIMP DM with an arbitrary abundance. As the abundance is inversely proportional to σv , an upper limit on the astrophysical flux can be used to set a lower limit on σv (see Eq. (19)). If WIMPs make up only part of the DM mass budget, its cross section is no longer restricted to be thermal. More generally, once the lower limit on the WIMP cross section exceeds the unitarity bound, fundamental WIMPs of this mass (DM or otherwise) will be totally ruled out [51].

VIII. TOWARDS CLOSING THE WIMP WINDOW
Future progress for decisively probing the WIMP paradigm requires improvement in indirect WIMP searches.
For Fermi, the most important thing for progress is the discovery of new dwarf spheroidal galaxies, especially any that are closer to Earth. The Dark Energy Survey (DES) is well poised to achieve this goal [143,144], that could lead to an improvement in limits on the single channels by about an order of magnitude [145]. New, more powerful gamma-ray instruments with better angular resolution or greater sensitivity are needed in the 10+ GeV range, such as GAMMA-400 [146,147] or HERD (High Energy cosmic-Radiation Detection) [148,149]. Sub-GeV probes such as PANGU (PAir-productioN Gammaray Unit) [150,151], AMEGO (All-sky Medium Energy Gamma-ray Observatory) [152], or ComPair (Compton-Pair Production Space Telescope) [153] will lead to better background subtraction for higher-energy searches. Otherwise, the sensitivity reach for Fermi will increase with the square root of exposure time. So far, Fermi has collected nearly 10 years of data, but has reported constraints on dSphs using 6 years of data [63].
AMS, on the other hand, sets bounds using a shorter exposure time of ∼ 2.5 years [64], and so can expect a greater increase based on exposure time alone. A key issue facing AMS limits is CR background/propagation uncertainties, which substantially restrict the strength of our conservative limits. Working towards a better theory understanding would significantly aid progress for decisively excluding the WIMP (see also Ref. [154]). Experimentally, better measurement of the local DM density will have a substantial impact (see Fig. 7), allowing a much stronger constraint on the lower DM mass. Analysis of new Gaia data will be vital here [155,156]. Future probes of the CMB only expect an improvement of a factor of a few. This is due to a fundamental bound of cosmic variance: we only have one universe to measure.
There are good prospects for improvements at high mass. The Cherenkov Telescope Array (CTA) is often considered to be decisive for higher WIMP masses (above ∼100 GeV). While it will have an important role in testing WIMPs, it will not be able to close the current WIMP window. Progress is required below these masses first, to make decisive statements about the status of the WIMP. Other than CTA [157], current generation telescopes, H.E.S.S. [86,87,89,158], VERITAS (Very Energetic Radiation Imaging Telescope Array System) [159,160], MAGIC (Major Atmospheric Gamma Imaging Cherenkov Telescopes) [161][162][163], HAWC (High-Altitude Water Cherenkov Observatory) [164][165][166][167][168], DAMPE (Dark Matter Particle Explorer mission) [169], and in future LHAASO (Large High Altitude Air Shower Observatory) [170,171], will also aid eventually closing up to the unitarity limit at ∼100 TeV.
As muon-dominated final states are consistently the least constrained, anything that improves constraints on muon-rich models is well motivated by our study to close the thermal WIMP window. A recent study of DM annihilation in the Milky Way Halo [107] did not include inverse Compton scattering (ICS), which is sensible when being conservative, because the Galactic interstellar radiation field and magnetic field are not well understood. However, as a consequence, the constraints presented on muon-heavy models are not very strong. A similar approach with a careful accounting for ICS could set stronger constraints on muon-rich models.

IX. CONCLUSIONS
Recently, there has been some growing sentiment that thermal WIMPs are on death row. However, such statements are often based upon direct detection scattering rates, or collider missing-momentum searches. In both cases, there is no well-defined and predictive scale for WIMP-SM interactions, and only specific aspects of any such interactions are being probed. Furthermore, any interpretation of such limits requires model-dependent choices and additional assumptions that cannot be model-independently related to the WIMP thermal relic cross section. Such limits may exclude model-dependent possibilities, but reveal nothing about whether the remaining possibilities are viable or not. Probes of WIMP annihilation products are tied to the fundamental nature of the WIMP as an annihilation relic, and so allow for a decisive statement about the exclusion of the WIMP.
We have calculated the first largely model-independent upper limit on the total WIMP annihilation cross section from data, meaningfully combining bounds from all kinematically possible DM annihilation products. We find that thermal WIMPs with s-wave 2 → 2 annihilation to visible final states have a lower exclusion bound of m χ 20 GeV. For the bound near 20 GeV, we have shown that the limit barely extends past the thermal relic scale, and pushing to higher statistical significance weakens this further. Enforcing a clearance factor of 2 leads to a lower limit on the DM mass of about 6 GeV, and a clearance factor of 10 gives a lower limit of about 2 GeV. For sub-GeV WIMP DM, the bound is severe, and extends to much lower masses.
The only way to decisively test the thermal WIMP scenario is to gain sensitivity to the total annihilation cross section down to the theoretical expectation. We have established an upper limit on the cross-section based on data, for the most generic WIMP scenario. This features an s-wave 2 → 2 annihilation cross section into visible final states, and a standard thermal history, with a radiation-dominated early universe. Other complications are possible, but they generally weaken limits. Considering a larger set of possibilities generally will just push the lower mass limit lower, which supports our point that GeV-scale WIMP DM is not even slightly dead.
We have discussed important improvements for moving towards closing the WIMP window. Discovery of new dwarf galaxies, with the aid of DES, would enhance limits from Fermi. This could significantly improve sensitivity, leading with photon-rich final states. Refinement of AMS results through better understanding of CR propagation uncertainties will provide a better probe of leptonic final states. Together, using our framework, these experiments or similar have potential to modelindependently exclude the generic WIMP up to the 100 GeV scale. Once this scale is reached, CTA will play a key role in excluding (or discovering) the WIMP. The remaining WIMP window is finite, and can ultimately be fully probed in a largely model-independent way.