Gamma-ray constraints on dark-matter annihilation to electroweak gauge and Higgs bosons

Dark-matter annihilation into electroweak gauge and Higgs bosons results in $\gamma$-ray emission. We use observational upper limits on the fluxes of both line and continuum $\gamma$-rays from the Milky Way Galactic Center and from Milky Way dwarf companion galaxies to set exclusion limits on allowed dark-matter masses. (Generally, Galactic Center $\gamma$-ray line search limits from the Fermi-LAT and the H.E.S.S. experiments are most restrictive.) Our limits apply under the following assumptions: a) the dark matter species is a cold thermal relic with present mass density equal to the measured dark-matter density of the universe; b) dark-matter annihilation to standard-model particles is described in the non-relativistic limit by a single effective operator ${\cal O} \propto J_{DM}\cdot J_{SM}$, where $J_{DM}$ is a standard-model singlet current consisting of dark-matter fields (Dirac fermions or complex scalars), and $J_{SM}$ is a standard-model singlet current consisting of electroweak gauge and Higgs bosons; and c) the dark-matter mass is in the range 5 GeV to 20 TeV. We consider, in turn, the 34 possible operators with mass dimension 8 or lower with non-zero s-wave annihilation channels satisfying the above assumptions. Our limits are presented in a large number of figures, one for each of the 34 possible operators; these limits can be grouped into 13 classes determined by the field content and structure of the operators. We also identify three classes of operators (coupling to the Higgs and $SU(2)_L$ gauge bosons) that can supply a 130 GeV line with the desired strength to fit the putative line signal in Fermi data, while saturating the relic density and satisfying all other indirect constraints we consider.


I. INTRODUCTION
Until the nature of dark matter and dark energy is understood, the remarkable success of the standard model of cosmology in accounting for observations will be less than completely satisfying. The most popular explanation for dark matter (DM) is that it is an early-universe relic in the form of an undiscovered, electrically neutral, nonhadronic, stable species of elementary particle. Of the many possibilities for the origin and nature of this purported new elementary particle, one that is particularly amenable to experimental and observational tests is that the new species is a cold thermal relic.
The cold thermal relic scenario assumes that the new particle species was established in local thermodynamic equilibrium when the temperature of the universe was larger than the mass of the particle. The processes establishing the initial equilibrium abundance are assumed to be the annihilation of the new species into standard-model (SM) particles and the production of the new species from initial SM-particle states. Then, as the universe cools to temperatures below the dark-matter mass, the relative abundance of the new species (relative to, say, the entropy density) "freezes out" as the rate of processes keeping the species in equilibrium falls below the expansion rate of the universe. The relative freeze-out abundance, and therefore the present mass density, is thus related to the cross section of dark-matter annihilation into SM particles.
Observationally the ratio of the present average dark-matter mass density to the critical density is well determined to be Ω DM h 2 0.12 [1]. The requirement that this mass density is saturated by the mass density of the cold thermal relic determines the annihilation cross section (and, in some cases, also the mass) of the new species. Since the mass and annihilation cross section required are both of order the weak scale, the cold thermal relic is known as a weakly interacting massive particle, or WIMP.
A key feature of the WIMP hypothesis is the requisite DM-SM coupling, and its close relationship to the present dark-matter mass density. Indeed, it is knowledge of the coupling of WIMPs to SM particles that provides the experimental and observational avenues that could lead to the discovery of the nature of the dark matter: detection More importantly, our results demonstrate that it is possible to set interesting limits on dark matter annihilations by combining various channels. They can also be used to infer additional predictions and learn about the nature of dark matter if a signal is observed in a particular channel.
We also identify an interesting effect for fermion DM coupling via pseudo-vector J DM where the p-wave annihilation cross section is enhanced (relative to the s-wave) at large M by coupling to the longitudinal mode of a massive vector boson in the final state. This actually has the effect, since saturation of Ωh 2 fixes the total cross section at freezeout, of suppressing the s-wave (i.e. present-day) annihilation cross section thereby alleviating constraints on these operators for large M .
Since there have been many works that discuss indirect-detection limits, it is appropriate for us to discuss why our analysis is new. Our starting point is the assumption that low-energy WIMP annihilation processes can be described by an effective field theory with the assumptions discussed above. Working within that framework, for each possible operator, the full SU (2) L × U (1) Y gauge invariance of J SM forces us to consider all possible final states. 2 Moreover, the requirement of gauge invariance determines the relative strength of possible annihilation channels. We expect the existence of multiple annihilation channels to offer much more discriminating power in the analysis. A related earlier work is Cotta et al. [22], where a more limited set of operators (without including the Higgs and without requiring full SU (2) L gauge invariance) was considered. In addition, some of these operators have been discussed in Refs. [23][24][25], with an emphasis on explaining the 130 GeV line as well as a possible additional line at around 114 GeV [26].
We emphasize that, in general, the EFT approach must be taken with some caution; the case at hand is no exception. Every relevant process involving WIMPs has a characteristic scale for the momentum transfer, and the accuracy of the EFT approach can be measured by the comparison between this scale and the suppression scale of the effective operators [27][28][29]. In this sense, the approach would become less accurate in describing the freeze-out and present-day DM annihilation if the suppression scale is close to the dark-matter mass. As it happens, in order to give the correct thermal relic abundances, some of the operators we study are forced to have Λ ∼ M and, in principle, one therefore needs to delve into possible UV completions in these cases. This caveat notwithstanding, given the large number of possible operators and final states we consider, we find it useful to use 'naïvely' the effective operator approach to give rough estimates. Considering the large systematic astrophysical uncertainties in the indirect detection measurements we utilize, a more accurate description within a UV complete model would not change our conclusions qualitatively.
The remainder of this paper is arranged as follows. In Sec. II we discuss some theoretical background on the WIMP relic density, DM halo profiles, annihilation photon fluxes, and continuum photon spectra from various DM annihilation channels. In Sec. III we discuss the individual experimental limits we have worked with, summarizing how they were obtained and how we have interpreted them. The large number of plots giving the results of our work are given in Figs. 4 to 37 in Sec. IV. We discuss the results and conclude in Sec. V. There are a number of appendices: Appendix A discusses in some detail the contribution to the photon spectrum due to inverse Compton scattering. Appendix B summarizes the magnitude of other systematic errors we did not take into account. Finally, Appendix C discusses some further analysis we have performed on one of the literature sources we have consulted for the case where there are two partially resolved lines in the photon spectrum.

A. Relic density
It is a well known phenomenon [30] that a massive particle species capable of annihilation into less massive species, and which is initially in thermodynamic equilibrium in the early universe, cannot maintain an equilibrium abundance to arbitrarily low temperature as the universe expands. Eventually, annihilations become too infrequent 3 owing to the dilution of the particles in the expanding volume of the universe, and creation of a pair of the new particle species becomes rare because it becomes exponentially rare to have collisions with sufficient center-of-mass energy. This leads to a "freeze out" of the abundance of the new particle species to some relic abundance. We shall assume that the measured average dark matter density of the universe, Ω DM h 2 = 0.12 [1], is entirely ascribable to the WIMP, which we will assume to be either a complex scalar or a Dirac fermion.
The phenomenon of thermal freeze-out is well described by the Boltzmann equation [30], and given the necessary initial data and annihilation cross sections as a function of temperature, one can simply solve this equation numerically to find the relic abundance. There does, however, exist an approximate method which has the benefit of expressing I. Milky-Way DM halo profile normalizations. We utilize a band of normalizations corresponding to the central value of ρ ± 2σ limits consistent with microlensing and galactic rotational velocities given in Iocco et al. [31]. For generalized NFW and Einasto profiles, we assume rs = 20 kpc, and r = 8 kpc, and take the ρ values from Iocco et al. [31]; for the Isothermal profile we assume rs = 5 kpc, and r = 8 kpc, and use representative values for ρ . the relic abundance in closed form. We shall assume that a non-relativistic approximation to the annihilation cross section can be written in the form σv N R = a + bv 2 , where a and b are constants, v 2 = 6T /M , and · · · indicates the thermal average. 4 The DM relic density can then be given to ca. 5% accuracy by [15,30] Ω DM h 2 = 1.04 × 10 9 GeV −1 M Pl g * (x F ) x F (a + 3b/x F ) × 1 for self-conjugate DM 2 for non-self-conjugate DM , (1a) (20) with T F the freeze-out temperature [30], M is the WIMP mass, M Pl = 1.22 × 10 19 GeV is the Planck mass, the numerical parameter c is chosen such that c(c + 2) = 1, and the number of relativistic degrees of freedom at freeze-out is taken to be g * (x F ) = 106 [15,30]. For real or complex scalar DM, g = 1, and for both Majorana and Dirac fermion DM, g = 2. For the results in this paper we will assume non-self-conjugate DM, since some of the operators under consideration vanish for self-conjugate DM. Working in an EFT framework and assuming that the DM annihilation proceeds through only one EFT operator O(x) of mass-dimension d, the term in the Lagrangian responsible for the annihilation can be written L(x) ⊃ Λ 4−d O(x). It then follows that σv is proportional to Λ 2(4−d) . Thus, for a given operator, Ω DM h 2 is a function of M and Λ. Since we assume the particle is a WIMP accounting for 100% of the observed DM in the universe, for a given mass M we shall numerically solve Eq. (1) for the EFT scale parameter Λ to satisfy the relation Ω DM h 2 = 0.12. This requirement implies that for each operator, for a given mass M there are no free parameters in the low-energy annihilation cross section.

B. DM Density Profiles
Although thermal freeze-out is predicated on the effective shutoff of annihilation in the early universe, this does not imply that WIMP annihilation remains negligible for the remainder of the evolution of the universe. Structure formation leads to large local densities of both baryonic and dark matter, and this obviously implies that the annihilation rate Γ ∼ nσv may again become large enough for detectable levels of annihilation to occur in regions such as the Milky Way Galactic Center (GC), or in the dark-matter-rich Milky Way dwarf companion galaxies.
Understanding the spatial distribution ρ(r) = M n(r) of the DM mass density in the relevant regions is important in the interpretation of any putative annihilation signal. The N -body simulation community is actively studying TABLE II. Summary of the main sources used to compute cross section limits, with corresponding choices of DM halo profile normalizations and J factors where applicable. [The J factor is defined in Eq. (6b).] If no value for a parameter was explicitly stated in a given reference, we indicate why we have assumed the indicated value. The region of interest (ROI) used in the analysis is indicated in parentheses following the notation of the source referenced. Further details in each case can be found in the referenced section number. In the final column we give the rescaling factor K applied in order to map published limits onto the halo profiles and normalizations given in Table I; further details on the rescaling are given in Sec. II C 1. For the most part the rescaling comes from scaling J factors with ρ 2 . The only exception is the isothermal case used by Weniger, where there is a significant difference in rs that affects the J factor by a factor of 2. a Not explicitly stated in the reference, but using this value correctly reproduces their J R3 , for which no point-source masking was applied. The reference also cites Catena and Ullio [33] where r = 8.2 kpc is mentioned. b The ROI for all cases is a 1 • radius circle centered on the GC, with |b| < 0.3 • excluded and no point-source masking applied. c J factors are not given explicitly in the reference, but we have verified these reproduce the limits. d Not explicitly stated in the reference, but assumed since this is the value used in Iocco et al. [31] from which the reference took all its normalizations.
this, but there is as yet no consensus in the literature for the exact form this so-called halo profile takes in galaxies such as our own. Consequently, it is conventional to quote all results assuming a number of different profiles. For constraints from Milky Way GC observations, we will utilize the (generalized) Navarro-Frenk-White, the Einasto, and the Isothermal profiles [16-18, 31, 32]. The generalized Navarro-Frenk-White (NFW) DM density profile is given by [31] ρ(r; γ) = ρ 0 (r/r s ) where the scale parameter is chosen to be r s = 20 kpc. The parameter γ controls the central slope of the profile: γ = 1 for the canonical NFW profile, while γ > 1 defines a profile with a steeper central region (where ρ ∼ r −γ ), known as a contracted profile (NFWc) which may arise due to gravitational interactions between the dark and baryonic matter as the latter cools during galaxy formation [34]. In either case, the profile is peaked toward r = 0 and "cuspy." The overall scale ρ 0 is chosen for galactic measurements such that the local (in the vicinity of the Sun) DM density ρ(r ) ≡ ρ takes some desired value (see Table I), while for extragalactic (dwarf) measurements, ρ 0 is chosen to obtain the correct observed stellar rotational velocities [19]. The Einasto profile is given by [31] where the scale parameter is chosen to be r s = 20 kpc. In line with large-scale numerical simulations, the parameter α = 0.17 is used. This profile is cuspy. Overall normalizations are chosen similarly to the NFW profile; see Table I. The Isothermal profile is given by [31] ρ(r) = ρ 0 1 + (r/r s ) 2 (4) where the scale-parameter r s = 5 kpc is chosen. This profile is known as "cored" since ρ(r) flattens off to a constant value at small r. Overall normalizations are again chosen similarly to the NFW profile; see Table I. We have scaled all the limits presented in this paper to reflect the halo profile parameters presented in Table I (original parameters assumed in the literature are shown in Table II; see Sec. II C 1 for more details on the rescaling).

C. Annihilation Photon Flux
For all the EFT operators under consideration, the detectable annihilation signal in γ-rays includes a) one or more monochromatic prompt photon "lines" arising from the direct annihilation DM DM → γX where X can be any SM boson consistent with gauge symmetries (X = γ, Z, h), and/or b) a prompt diffuse continuum photon flux arising from final-state radiation or radiative hadron decay arising when considering various non-photon primary DM annihilation products: DM DM → XY → γ's where X and/or Y is not a photon.
For each particular operator there will be a characteristic per-annihilation differential prompt photon spectrum where BR f ≡ σ (DM DM → f ) / σ total is the branching ratio for the annihilation mode to final-state f , dN f γ (E)/dE is the per-annihilation differential prompt photon spectrum for that mode, and the sum runs over all annihilation modes for the operator. Given this spectrum, the differential DM-annihilation prompt photon flux which an experiment would observe at photon energy E is given by [32] dΦ Hereafter we will assume non-self-conjugate DM, i.e., complex scalars and Dirac fermions. In the equation above, ρ(r) is the DM density profile at Galactocentric distance r. The astrophysical "J factor" [Eq. (6b)] is computed by integrating the squared DM density along the line of sight (LOS) from the Earth, and over the relevant Galactic lat/long coordinates. 5 The integral over b and l defines a search region of interest (ROI). A point observed at LOS distance s at lat/long coordinates (l, b) has Galactocentric distance r = s 2 + r 2 − 2sr cos l cos b 1/2 where r is the Sun-GC distance. We will always take r = 8 kpc [31]. Finally note that the polar angle θ from the GC satisfies cos θ = cos l cos b.
It is our aim to place upper limits on the non-relativistic cross sections σv using experimentally determined upper limits on such photon fluxes.

Profile rescaling
The constraints derived in this paper draw on published gamma-ray line searches and limits on diffuse emission from several sources, each of which employ different assumptions for Milky Way halo profiles. We give the parameters used by these references in Table II. For consistency, however, we would like to derive limits for a single set of halo profile parameters (given in Table I).
The difference in halo profiles is encapsulated in the J factor; typically we wish to take some given J factor from the literature (using some values of r , ρ and/or r s ) and rescale it to values of r , ρ and/or r s that are our choices. All other things being equal, the limit for σv ∝ 1/J, so the cross section limit that would be set with the new set of halo parameters is Here J is the J factor computed with the set of parameters we wish to use and J is the J factor we compute using the parameters in the literature. We thus rescale the published limits by the factor K (see Table II). We find that scaling the J factors with ρ 2 accounts for most of the difference in the halo profile normalizations. This is primarily because the variation in ρ values used in the literature is large (up to 40%), while the variation in the r values is less than 10%, and the r s values are almost uniform. There is only one case where there is a large difference in r s ; for the isothermal profile for Reg4 of Ref. [20] where r s = 3.5 kpc → 5.0 kpc. In this case the rescaling factor was roughly a factor of 2 smaller than one would have obtained from simply scaling J with ρ 2 , accounting for the fact that the core is less dense when r s is increased.
Note that for the Fermi ROIs, we have computed these rescaling factors without considering point source masking effects which are estimated to impact the J factor by less than about 10%. Our justification for this simplification is that the halo profiles are smooth away from r = 0, so that the ratios J /J computed with or without point-source masking being considered should be very nearly equal (certainly much closer than 10%).
We also point out that the rescaling of the Reg4 values from Ref. [20] is appropriate even though the ROIs in that reference were chosen based on an expected signal/background analysis with a specific profile in mind. Weniger [20] reports cross section results from each such ROI choice assuming a number of different profiles in the conversion from flux to cross section while holding the ROI fixed, and it is only these flux-to-cross-section conversions we are rescaling. Similarly, the ROIs in Ref. [16] were chosen based on an expected signal/background analysis; but again, we are simply holding their chosen ROI, and hence fluxes, fixed and are only rescaling the flux-to-cross-section conversion.

D. Photon Spectra
In the models we study, dark matter annihilation produces primary electroweak gauge bosons (γ, Z, W ± ), higgs bosons, or fermions. In addition to primary photons, the subsequent shower and hadronization of final states will produce abundant hadrons, giving rise to additional prompt photons dominantly through decay of neutral pions π 0 → γγ. To extract cross section upper limits from the measured flux upper limits requires knowledge of this perannihilation prompt photon spectrum, see Eq. (6a). Here we summarize the physics of prompt photon production from DM annihilation.
We begin our discussion with monochromatic photons. The γγ and γh final states give rise to photon lines. 6 The contribution to the photon spectrum is a delta function: for a γX final state, and N γ counts the number of monochromatic photons in the decay mode (N γ = 2 for the γγ annihilation modes and N γ = 1 otherwise). The photon line from a γZ annihilation mode is similarly peaked at E γ = M − m 2 Z /4M , but obviously has a finite width due to the finite width of the Z. The width of the line will be relevant to our discussion of experimental constraints later and is shown in Fig. 1.
For annihilation modes W + W − , ZZ, Zh, hh, ff , γh, and γZ, there is also a continuum prompt photon spectrum. We employ Pythia 8.176 [36] to perform the shower and hadronization (and hadron decay) to obtain the photon spectra. Following Ref. [37], we create a fictitious spin-1 resonance at m = 2M with a width of 1 GeV, which we populate using fictitious e + e − colliding beams free of initial-state radiation at √ s = 2M . We set by hand the resonance decay channels to be, in turn, 100% to each of the desired primary final states, and the per-annihilation differential photon spectra are computed by correctly normalizing 7 the histogrammed photon energies from 10 5 simulated annihilations per mass point per primary final state. Good agreement is found between these computed Pythia spectra and a set of spectra from Cirelli et al. [37] for DM masses below 300 GeV. For larger DM masses, the effects of electroweak (EW) final-state radiation (FSR) of W s and Zs (which are not included in the Pythia shower) can significantly alter the spectra at low E γ [37,38]; these effects have been included in the spectra in Cirelli et al. [37]. However, we shall see below that we will be primarily interested in the spectrum for E γ ∈ [10, 100] GeV at large DM mass. In this photon energy range, the EW corrections are unimportant for M 1 TeV except for primary final states containing charged leptons, which in any event give rise to about an order of magnitude fewer photons per annihilation compared to other channels. Furthermore, the branching fraction for annihilation to leptonic final states in particular, and fermionic final states more generally, is subdominant at large M (see the result plots in Sec. IV), with the exception of operator XX-1 where all annihilations are mediated through an s-channel Z or γ. Therefore, since the error of neglecting the EW FSR is small, and since the necessary spectra are not all available in Cirelli et al. [37], we have chosen for the sake of uniformity to use the Pythia spectra which we have computed rather than attempting to use the corrected spectra from Cirelli et al. [37].
1. Qualitative understanding of the shape of the photon spectra The shape, and in particular the DM-mass dependence, of the continuum photon spectrum differs qualitatively depending on whether the primary final states are leptonic, non-leptonic and color neutral, or colored. We consider illustrative cases for the latter two: annihilation to W + W − and annihilation tobb (light qq modes are qualitatively similar), and comment on the leptonic modes.
Consider first the W W case. In Fig. 2 we plot the photon spectrum decomposed into photons from π 0 decay and those from final state radiation (FSR) of charged leptons (omitting other hadronic-decay photons for clarity). For E γ 1 MeV, the production mechanism is dominantly due to pion decay. The strongest continuum γ-ray constraints are for energies from 0.1 GeV to 100 GeV, so we will always be in the regime where the behavior of the photon spectrum from diboson annihilation channels is dominated by photons from π 0 → γγ. Meanwhile, the FSR from leptons has a spectrum that goes like dN γ /dE γ ∼ E −1 γ and is subdominant for E γ 1 MeV. Since W is color neutral, the spectrum and multiplicity of photons is determined by the independent decay and shower of each W , boosted into the center of mass frame of the annihilation. As M increases, the W s are increasingly boosted and the decay products also become more energetic while the multiplicity (namely the pion yield) remains roughly constant. We can observe this in the shifting of the photon spectrum in Fig. 2 to higher energies. Furthermore, this migration to higher energy decay products, while the yield remains constant, leads to the effect that as M increases there is a decrease in the number of photons produced at low energy. This can also be observed in e.g., Panel B of The photon spectrum arising from DM annihilation, decomposed according to the production mechanism of either π 0 decay or final state radiation of charged leptons. We show annihilation to W W (Panel A) and to bb (Panel B), for DM mass of 130 GeV and 800 GeV. It is clear that for Eγ 1 MeV, π 0 decay dominates the spectrum. In the case of the W W final state, the shift of the pions to higher energies reflects the increasing boost of each W (and thus of its decay products). For the bb final state, the yield of pions increases over the entire energy range of interest. We have omitted the spectra from other hadronic decay modes for clarity.
around E γ ∼ 0.1 GeV as M increases.
Since most of the operators considered here have annihilation modes into gauge bosons or higgs, the mass dependence of the γ-ray spectra has similar behavior to the W W case, modulo the effect of changing branching fractions. Some operators which have only diboson annihilation modes also show first an increase in the low-energy photon spectrum with increasing M below or around a hundred GeV before the decrease discussed above sets in. This has less to do with the intrinsic photon spectrum from each annihilation mode and more to do with the branching fractions for each mode. If the γγ mode is present, then the branching fraction to γZ starts off small around M ∼ m Z /2 and rapidly increases with M giving the observed increase in the continuum spectrum (e.g., Operator VI-1 in Fig. 4 of Sec. IV). In this case, there may be a W W or ZZ mode entering once M gets to m W,Z , which results in a further increase since the spectra are normalized per-annihilation and twice as many W/Zs per annihilation gives twice as many photons per annihilation (the W and Z intrinsic photon spectra are almost identical, modulo kinematic effects). For M above a (few) hundred GeV, the low-energy spectrum begins to decrease again as per the discussion in the previous paragraph. Explanations of a similar nature are possible for cases where other mixtures of diboson annihilation modes are present.
For annihilation to bb, shown in Fig. 2, the dominant photon-production mechanism is still neutral pion decay, but the same qualitative picture does not obtain: there is a significant increase in the pion yield with increasing M because in this case the initial hard scale for the parton shower is set by M and so the number of FSR gluons and other strongly interacting partons produced in the shower increases significantly as the DM mass increases. Therefore, although the average pion energy does increase, the increase in the pion yield at the same time means that the pion yield at low energies actually slowly increases as M increases; the photon yield thus shows a similar increasing trend. This is evident in the slow overall monotonic increase in the spectrum as M increases for Operator XX-1 shown in Fig. 37. (For the other operators considered here, this behavior gets cut off near M ∼ 100 GeV since the branching fraction to quarks typically becomes subdominant once diboson annihilation modes become kinematically allowed.) Finally, for annihilation to charged leptons, the dominant photon production mechanism is direct final-state photon radiation from the leptons. This spectrum falls roughly as E −1 γ over most of the energy range with a high-energy kinematic cutoff for the photon spectrum at M , while the normalization of the photon spectrum grows slowly as M increases. Again however, with the exception of operator XX-1 already noted above, fermionic final states tend to have rapidly falling branching fractions once on-shell diboson final states become kinematically allowed as M increases, so the increase in the spectrum at low E γ ∼ 0.1 GeV is cut off around M ∼ 100 GeV.

III. LIMITS
In this section we discuss the experimental limits used in constraining the effective operators. We consider line searches in gamma-ray data from Fermi-LAT and H.E.S.S., as well as constraints on diffuse gamma-ray emission derived from Fermi-LAT observations of the GC and dwarf galaxies, for which the continuum photon emission is relevant. We describe our methods of extracting limits for the operators, re-interpreting published results for the specific final states and branching ratios of each operator. Our final results are shown in the figures in Sec. IV, along with a universal caption and summary of the constraints on the operators.

A. Low-energy line limits: Fermi-LAT
The Fermi-LAT Collaboration [16] presents 95% confidence level (CL) upper limits (UL) on the flux Φ of photons from regions on the sky centered on the Galactic Center due to a strictly monochromatic underlying photon line in the photon energy range 5 to 300 GeV. The underlying photon energy E γ was scanned across this range, and for each energy they obtained a flux limit by performing a maximum likelihood analysis in which the Fermi-LAT line shape plus a power-law background was fitted 8 to their 3.7-year data over sliding energy ranges E γ ± 6σ E where σ E is the Fermi-LAT energy resolution (68% containment half-width). The flux limit Φ they report is the 95% CL UL of the normalization of the line in this fit, multiplied by a factor giving the effective exposure time of the experiment.
The Fermi-LAT search methodology selected a different signal region of interest (ROI) for each of four DM halo profiles (NFW, Einasto, Isothermal, NFWc γ=1.3 ) they considered in order to maximize the expected signal-to-background ratio 9 as estimated through Monte Carlo simulations. Thus, they present four distinct flux limits, one for each ROI. These ROIs are denoted Rα where α is the angular size of a circle centered on the GC, with the region |l| > 6 • and |b| < 5 • masked; regions around known γ-ray point sources are also masked for all ROIs except R3 [16].
We derive constraints on σv from the flux limits presented in Ref. [16], which are given at discrete E γ in the energy range 5 to 300 GeV; we have interpolated (linearly in log Φ vs. E γ ) these limits to intermediate energies where necessary. If the underlying photon line arising from a γX final-state mode is monochromatic or monochromatic to an excellent approximation, we could merely use the flux upper limit at E γ = M − M 2 X /4M combined with Eq. (6b) and the J factors in Ref. [16] (rescaled per Table II) to derive limits on σv: where N γ is the number of photons in the line per annihilation: 2 for γγ and 1 otherwise. This is the case for the γh annihilation mode; complications however arise for the γZ and γγ modes for the operators we consider. The photon line arising from γZ annihilation modes has an intrinsic width of order 1 GeV, which can be significant with respect to the ca. 5-10% Fermi-LAT energy resolution at low energies [16,39] (see Fig. 1). However, even for a line with intrinsic width of 50% of the Fermi-LAT energy resolution, the resulting observed spectral feature is only expected to be broadened by around 11% after convolution with the Fermi-LAT energy dispersion function. If such a feature were fitted with a line shape derived from the assumption of a purely monochromatic intrinsic line, the number of photons would be underestimated by about 7% (Ref. [16], App. D2-3) which would set the flux limit correspondingly more stringent than it should. However, this effect is much smaller than both the expected statistical fluctuation (ca. 50%) in the limits set by Ref. [16] and the uncertainty in the halo profile normalization (a factor of ca. 10; see results) and it is roughly on the same order as other systematic effects (ca. 10%) which are present in the Fermi-LAT analysis and ignored in their flux upper limits [16]. Provided the linewidth is expected to be less than about half the energy resolution, we therefore make no correction for the finite Z width, and merely use the Fermi-LAT flux limits derived under the assumption of a monochromatic intrinsic line in deriving limits on [σv] γZ . To satisfy this linewidth requirement, we follow the older Fermi-LAT line search analysis [32] and only present limits on σv from the γZ annihilation mode for E γ ≥ 30 GeV, which corresponds to M 63 GeV (see also Fig. 1).
There is a significant complication [40] in interpreting the flux limits for all the operators with γγ annihilation modes: due to their gauge structure, all such operators under consideration also have a γZ annihilation mode which yields a second line at an energy δE = m 2 Z /4M lower than the γγ line. Since the unbinned likelihood analysis in Ref. [16] is sensitive to both shape and normalization of the line spectral feature, this makes the interpretation of the single-line flux limits difficult in the region of DM masses where the two photon lines are separated by an amount of the order of the Fermi-LAT energy resolution. Depending on the WIMP mass, we have adopted a few different techniques to interpret the single-line flux limits for these operators. In all cases, we choose to use the flux limits to set a single limit on the quantity (2σ γγ + σ γZ ) v.

M < 63 GeV
This mass range corresponds to E γ ≤ 30 GeV for the γZ channel (when kinematically allowed), so we do not set a γZ line limit in this mass range. We find the 95% CL UL on σ γγ v using Eq. (8) with the flux at E γ = M , and use this to set, at each WIMP mass, the limit where the cross section ratio in the [ · · · ] brackets is computed using the analytic expressions in CKW and rescales the σ γγ v limit to the quantity (2σ γγ + σ γZ ) v.

63 GeV < M < 80 GeV
In this mass region, we have two lines which are relatively well separated 10 assuming a 5 to 10% energy resolution for Fermi-LAT [39]. We thus use Eq. (8) to set independent 95% CL UL on σ γγ v using the flux at E γ = M , and on σ γZ v using the flux at E γ = M − m 2 Z /4M . We use these two quantities to set, at each WIMP mass, the limit where the cross section ratios in the [ · · · ] brackets are computed using the analytic expressions in CKW.

80 GeV < M < 160 GeV
In this intermediate mass region, the spectral feature expected in the Fermi-LAT data corresponding to the two lines will be very broad with respect to the energy resolution, and may be a double-peaked structure depending on the relative strengths of the lines. Such a feature looks very different from the line shape used in the Fermi-LAT line search [16] and since, as noted above, that analysis was sensitive to the shape of the fitted feature as well as its normalization, we exercise an abundance of caution and do not set any limits on σv in this region. This caution notwithstanding, we do not expect the limits to differ from those at nearby energies by more than a factor of a few. Since a reproduction of the analysis in Ref. [16] with a different assumed line shape is a task best suited for the experimental collaboration, we would encourage the Fermi-LAT Collaboration to look into setting flux limits for the case where there may be two partially resolved lines present, with some variable ratio of strengths.

160 GeV < M
In this mass region, the two lines are sufficiently close that the spectral feature expected in the Fermi-LAT is merely a slightly broadened (width increased by 10%) line and a single limit can be set on the combined quantity (2σ γγ + σ γZ ) v. We do this by using the flux upper limit from Ref. [16] evaluated at the flux-weighted average energȳ where the cross sections here are computed using the analytic expressions in CKW. The limit is set as B. High energy line limits: H.E.S.S.
The H.E.S.S. Collaboration [17] presents 95% confidence level (CL) upper limits (UL) on the photon flux Φ from a strictly monochromatic photon line in the energy range 500 GeV to 20 TeV. The ROI is a 1 • radius circle centered on the GC, with |b| < 0.3 • excluded. No masking of point sources was performed. Their flux upper limits were extracted by performing a maximum likelihood analysis in which a Gaussian peak (photon line) was fitted along with a parametrized background to 112 hours of data from the H.E.S.S. VHE γ-ray experiment. The position of the Gaussian peak was scanned over the search energy range (with the standard deviation constrained to be equal to the H.E.S.S. energy resolution; see below), and the Gaussian normalization and background parameters were fitted at each search energy. The line flux limits reported are the 95% CL UL on the fitted Gaussian normalizations.
Flux limits are only presented at discrete E γ on the energy range 500 GeV to 20 TeV in Ref. [17]; we have interpolated (linearly in log Φ vs. log E γ ) these limits to intermediate energies where necessary. Isothermal limits here are very weak as a consequence of the cored profile and the ROI which is restricted to a very small region near the GC. It was necessary for interpreting the H.E.S.S. limits to compute the J factors defined in Eq. (6b) for the halo profiles of our choice, given in Table I. We did this for by computing J factors for the parameters used in the publication (see Table II) and then rescaling the resulting J factors as necessary, as discussed in Sec. II C 1.
Additionally, in computing cross section limits from the fluxes, we neglected the small shift in the photon line position for final states involving a massive particle (the largest this shift gets is about 7 GeV for the γh line at M = 500 GeV) and have simply evaluated all limits assuming E γ = M . Furthermore, since the H.E.S.S. energy resolution is 17% at 500 GeV, dropping to 11% at 10 TeV [17], it is always a good approximation to combine the unresolved γγ and γZ lines into a single spectral feature to derive combined limits on (2σ γγ + σ γZ )v for all the operators with both γγ and γZ annihilation modes. The width of the photon line from γZ has no significant effect here.
To summarize, we derived cross section limits from the H.E.S.S. flux upper limits using, as appropriate, C. Inclusive spectrum limits: Hooper et al. analysis Hooper et al. [18] presents 95% CL UL on the inclusive photon spectrum from the Galactic Center region based on the 3.7-year Fermi-LAT data in a model-independent form as 95% CL UL limits on the quantity F bin j ≡ (σv/M 2 ) bin j dN total γ /dE γ dE γ for the four photon energy bins of 0.1 to 1 GeV, 1 to 3 GeV, 3 to 10 GeV and 10 to 100 GeV. These limits were derived using a signal-template-based subtraction methodology with a photon-flux template proportional to LOS ρ 2 (r[s, l, b]) ds, rather than the ROI-integrated approach used in the other references we have used.
We employed these model-independent limits to set limits on [σv] total , taking the most stringent of the limits from the four bins listed above. For each M , the prompt photon spectrum dN γ /dE γ (including all annihilation modes) was numerically integrated over each of the four energy bins defined above. We have included the photons arising from lines where present. For each bin, the limit is formed (the factor of 2 is a re-interpretation of a limit set for self-conjugate DM in [18] For large M , bin 4 (10 to 100 GeV) typically gives the most stringent limit (cf. the comments in Sec. II D on the regions of the continuum photon spectrum in which we are interested). Note that the inclusion of line photons makes the resulting limits much more aggressive, compared to the case if only continuum photons were present, for M 100 GeV (the exact cutoff here depends on the identity of X for a γX annihilation mode). We have made no attempt to convolve any photon lines here with the Fermi-LAT experimental resolution; doing so would smooth some of the sharp transitions evident in our result plots at values of M where a line crosses one of the bin edges (see Sec. IV), but would make no other qualitative changes.
Ref. [18] also presents limits in a model-dependent form assuming 100% branching fraction to, in turn, W W , ZZ, bb, cc, µ + µ − , and τ + τ − final states. We find that there is a multiplicative factor of ca. 2.8 discrepancy between the model-independent (more stringent) and the more model-dependent (more conservative) limits once the continuum spectra are used to convert the model-independent limits to limits on the specific final states listed above. Following Hooper (private communication), we resolve this discrepancy in favor of the more conservative limits. We do this by multiplying all σv limits set for our operators using the model-independent F, as described in the paragraph above, by a factor of 2.8.

D. Continuum limits: Fermi-LAT
Assuming particular annihilation channels, the Fermi-LAT Collaboration [19] also presents 95% CL UL on the inclusive photon spectrum using observations of 10 Milky Way dwarf companions (the "stacked dwarf" limits). Their limits were obtained by performing a simultaneous maximum likelihood analysis on all 10 dwarfs in which the expected DM annihilation signal was fitted along with a diffuse galactic background model to 2 years of Fermi-LAT data for photons in the energy range 0.2 to 100 GeV. The normalizations of the two components were floated in the fit and cross section limits were extracted from the normalization of the DM signal component assuming in turn, that the annihilation was 100% via each of the modes bb, W W (equivalent to ZZ for these purposes), µ + µ − , and τ + τ − .
Unfortunately, since this analysis is sensitive to spectral shape, and since the continuum photon spectrum for annihilation modes have different shapes in general, it is not possible to simply reverse-engineer these limits to set the exact 95% CL UL considering all the continuum contributions from a combination of annihilation modes (which is the case for the operators we consider). We therefore extract approximate and conservative limits on [σv] total by requiring that the individual W W/ZZ, bb, µ + µ − , and τ + τ − contributions to σv for any operator do not violate the individual experimental upper bounds. For a given final state, we take the limit where [σv] 95% CL UL f is the limit taken directly from Ref. [19] and BR f is the branching fraction for the annihilation mode f ∈ {W W/ZZ , bb, µ + µ − , τ + τ − } computed using the analytic expressions in CKW. The factor of 2 is a re-interpretation of a limit set for self-conjugate DM in Ref. [19] to our choice of fermion or complex scalar DM.
These constraints from dwarf galaxies as we have implemented them are conservative and approximate and are not necessarily indicative of the full ability of such measurements to constrain these operators. Based on Fig. 14 of Ref. [18], we expect the exclusion reach for an analysis which took into account the entire continuum spectrum at once rather than only looking independently at each component of the spectrum from each final state and comparing to the relevant final state limit from Ref. [19] (thereby paying a cost of BR f ) would be similar to the inclusive galactic center limits from Ref. [18] in cases where there are no lines.
Finally, we note that as this paper was being finalized, the Fermi-LAT collaboration released updated constraints stacking 25 dwarf satellites based on 4 years of data [41]. These limits are weaker than expected, and in particular weaker than the limits we take from Ref. [19]. However, since these limits are not constraining (as implemented), the updated analysis in Ref. [41] ultimately does not affect our results. Cohen et al. [21] finds 5.5σ (local) evidence for a photon line in the 3.7-year Fermi-LAT data at roughly 130 GeV, and based on the same data, presents limits on the ratio R of the annihilation cross section via modes which give rise to continuum photons, to the annihilation cross section giving lines (γγ and/or γZ). These limits are based on observations of photon fluxes from Fermi-LAT in the energy range 5 to 200 GeV in an annulus around the GC with inner and outer radii 1 • and 3 • , respectively. We compare the two classes of R limits in Ref. [21] -shape and supersaturation -to the ratio R th = (σ W W v + σ ZZ v + σ bb v) / (2σ γγ v + σ γZ v), computed using the analytical expressions in CKW.
The supersaturation limits on R are derived without assuming a background model, requiring only that the photons from DM annihilation do not supersaturate the observed photon spectrum. These limits are only presented for the case where 100% of the continuum photon flux comes from annihilations going to a W + W − primary final state, but due to the similarity of the spectra, these limits are also applicable if the annihilation is to ZZ or bb. Although the limits are very conservative, they are more robust to details of astrophysical backgrounds or the DM model than the shape constraints.
We briefly summarize the procedure they followed to obtain the supersaturation limits: for each M , the number of line photons was found by performing a maximum likelihood analysis in which the photon lines for γγ and γZ were fitted to binned photon-count data, marginalizing over the relative normalization of the lines. The sum of the signal from the two lines was taken as the number of line photons. To find the upper limit on the continuum normalization, the energy bin where the continuum photon spectrum is expected to peak relative to an assumed power-law astrophysical background with a spectral index of 2.8 was first found (this is the only point where spectral information was used); for a pure W W annihilation mode, this was determined to be the bin 10 to 20 GeV. The normalization for the continuum annihilation was then scaled up until the number of photons in the 10 to 20 GeV bin violated the 95% CL UL from data. Correcting for the effective area in Fermi-LAT, this results in a limit on the continuum-to-line cross section ratio.
The second (much stronger) class of limits presented in Ref. [21] are shape limits, derived from a maximum likelihood analysis including the two photon lines, the continuum photons, and a single power-law background over the energy range 5 to 200 GeV. At each M , the fit was performed for the ratio of continuum-to-line photons in that energy bin, marginalizing over all other parameters including the relative strengths of the two lines. This procedure effectively resulted in a 2-dimensional likelihood profile in the M − R plane; the 2σ limit line in this M − R plane is taken as the 95% CL UL on R. We consider their results for each of the continuum-photon-producing final states in turn, 100% W W/ZZ and bb; the µ + µ − and τ + τ − modes were also considered there, but for our operators annihilation to charged leptons is typically subdominant.
It should be borne in mind that the analysis of Ref. [21] marginalized over the σ γγ v/σ γZ v ratio, while the operators here have fixed line ratio for a given M . Furthermore, the limits are most applicable if there is some region of M where a good fit to the photon line is possible; note that the line fits are discussed further in Sec. III G and also indicated in panel C of our results figures. With more recent data [16], evidence for the line has weakened, which requires fewer line photons and would lead to correspondingly weaker limits on the continuum-to-line ratio if this analysis were to be updated. Furthermore, Ref. [21] assumes 100% annihilation to each of the various continuum-photon-producing modes, while the operators we consider invariably have some combination of multiple continuum-photon-producing annihilation modes. Note also that no inverse Compton scattering (ICS) contribution was included in the continuum spectrum, as it should be sub-dominant for W W/ZZ and bb, and the contribution to the continuum spectrum of the Z in the γZ channel is not included. Finally, none of the limits presented is applicable to operators with γh final states as this annihilation mode introduces a new line which would dramatically improve the goodness-of-fit in the M ∼ 155 GeV region.
Nevertheless, even given all these issues, the limits remain indicative for the case where a least one of the annihilation modes γγ or γZ is present, and the operator's continuum photon spectrum is dominated by one of W W/ZZ or bb. Overall, one should interpret the limits qualitatively, not quantitatively: for operators with continuum-to-photon ratios near the limit, one would really have to be careful in re-doing the full likelihood fit with the correct known branching ratios for the various modes for the specific operator in question to definitively settle the issue of whether or not the operator is excluded.

F. Line-like feature: Weniger-like analysis
It is interesting to examine whether it is possible for the photon lines arising from these operators to explain, with the correct normalization, the line-like feature near 130 GeV 11 first reported at 4.6σ local significance by Weniger [20] in the 3.7-year Fermi-LAT data, assuming for these purposes that we take this feature seriously as a DM signal.
Assuming the annihilation is to γγ, Weniger [20] reports values for σv required to account for this feature. For a γγ mode, the DM mass must clearly be 130 GeV to explain the line. For the operators we consider here, there may be annihilation modes to any of the final states γX where X = γ, Z, h, where the DM mass must be M = 130, 144, and 155 GeV, respectively, to explain the observed line. Eq. (6a) indicates that we must account for this change in mass by performing a rescaling of the cross section result from Ref. [20]. To be explicit, where an operator has annihilation modes to γX (X = Z, h), we rescale the γγ result [σv] Weniger from Ref. [20] as follows: where one factor of 2 is due to the assumption of γγ and the other factor of 2 is a re-interpretation of a limit set for self-conjugate DM in [20] to our choice of fermion or complex scalar DM [see Eq. (6a)]. In addition, we perform the rescaling necessary to account for differing halo profile normalizations, as discussed in Sec. II C 1; the rescaling factors we use are given in Table II.
Owing once again to the only partial resolution in the Fermi-LAT data of the photon lines, we cannot naïvely use the required line cross section from Ref. [20] to give individual required cross sections for the γγ and γZ modes when both are present (see discussion in Sec. III A). Assuming either M = 130 GeV or M = 144 GeV, the secondary line at either 114 or 144 GeV, respectively, would give a significant photon contribution at the position of the line at 130 GeV and/or vice versa depending on the relative normalizations of the lines. Similarly, since the line separation here is still on the order of the Fermi-LAT energy resolution (even if assumed to be 10%), one also cannot assume the lines are completely unresolved to derive a single required combined value for (2σ γγ + σ γZ )v.
Nevertheless, while we acknowledge this issue, for the purpose of giving a rough estimate of the annihilation cross section which would be required for our operators with both γγ and γZ annihilation modes to explain the line reported at 130 GeV in Ref. [20], we have performed the following approximate analysis. We note that it is usually the case that one or other of the lines from either the γγ or γZ mode is dominant for M in the approximate range 120 to 150 GeV (in terms of photon flux, for operators VI-1, IV-2, VIII-1, or VIII-2, it is the γγ mode which dominates by about a factor of 4, while for operators VI-3, VI-4, VIII-3, or VIII-4, it is the γZ mode which dominates by about a factor of 2.5). Therefore, we shall make the approximation that the line at 130 GeV reported in Ref. [20] is entirely due to the dominant of the γX (X = γ, Z) modes. As noted above, there will be a second line in either case, but we ignore this complication in our approximate treatment here (see the next section for a slightly more careful treatment of the same issue making use of results from a different reference). We give results for the quantity (2σ γγ + σ γZ ) v, requiring that the normalization of the dominant mode can explain the line at 130 GeV. To be explicit, if the γX (X = γ, Z) annihilation mode explains the line at 130 GeV, we extend Eq. (17) above to the two-line case: The factors in the brackets [ · · · ] are computed using the analytical expressions in CKW and merely rescale the individual required σ γX v (X = γ, Z) cross sections for the line at 130 GeV to the quantity (2σ γγ + σ γZ ) v.

G. Line-like features: Cohen et al. analysis
In addition to presenting limits for the ratio of continuum-to-line cross sections, Cohen et al. [21] also presents evidence for the existence for a line in the Fermi-LAT 3.7-year data in an annulus centered on the GC with inner and outer radii of 1 • and 3 • , respectively, with local significance of 5.5σ relative to the null hypothesis of no line. To reach this conclusion, they performed a likelihood analysis fitting binned photon count data in the energy range 5 to 200 GeV to a single power-law background plus lines at E γ = M (with normalization N γγ ) and at E γ = M − m 2 Z /4M (with normalization N γZ ), both suitably convolved with the Fermi-LAT energy dispersion. No continuum photon contribution from DM annihilation was included in this analysis. In the fit, M and θ γZ/γγ = arctan (N γZ /N γγ ) were scanned over while the background normalization and spectral index, as well as the total line normalization N γγ +N γZ , were marginalized over. The best fit point for this analysis was found to be M = 130 GeV and θ γZ/γγ = 0. (An earlier analysis [26] fitting the data with γγ and γZ lines found very similar results; we have followed Cohen et al. [21] since that work also gave limits on the continuum-to-line ratio.) Cohen et al. [21] presents 1-, 2-, and 3-σ delta-log-likelihood contours for their fitting procedure in the M − θ γZ/γγ plane, which we have reproduced in Fig. 3, along with the theoretical values for θ γZ/γγ = arctan (σ γZ v/2σ γγ v) for the eight operators that have both γγ and γZ lines, as computed using the analytical expressions of CKW. We see that the relative normalizations of the lines for these operators are such there are regions around M = 130 GeV or M = 144 GeV where the sum of the two lines could fit the Fermi-LAT data very well (within 1σ in delta-log-likelihood). Cohen et al. also gives contours displaying the total normalization (photon count) N total ≡ N γγ + N γZ of the lines in the M − θ γZ/γγ plane which, knowing the ROI-averaged mission-time-integrated exposure-times-effective-area E ROI relevant to the Fermi-LAT data they considered, would allow us to extract the cross sections required for the operators with both γγ and γZ modes to explain the fitted lines. As E ROI was not relevant to the analysis presented in Ref. [21], it was not given there. However, we have made use of the Fermi ScienceTools package to reproduce a sufficient portion of the data-extraction and analysis in that reference to enable us to compute it ourselves (see Appendix C for the full details of our analysis and the cross-checks we have performed). We find that the ROI-averaged exposure factor for photon energies from ca. 120 to 150 GeV (applicable to the line-fits presented in Ref. [21]) is E ROI = 1.05 × 10 11 cm 2 s. Knowing this factor, and extracting the fitted total line-normalization N total from the contours in Fig. 2 of Cohen et al. [21], we have computed the photon flux corresponding to the fitted line normalization as Finally, making use of Eq. (6a) we have converted this to the annihilation cross section required as a function of M to explain the line for our particular operators having both γγ and γZ modes, for the (operator-specific) range of M within the 3σ |delta-log-likelihood| contours shown in Fig. 3. In performing this conversion, we have also computed the J factors for the ROI in Ref. [21], which we give in Appendix C.

IV. RESULTS
In Figs. 4 through 37 we present the indirect detection limits for all operators with s-wave annihilations as studied in CKW. The cross sections for various operators are listed in Tables VI-XX of CKW. We refer to a process studied in CKW by a table number, and an operator number. For instance, operator VI-3 is the third operator of Table VI in CKW, in this case, φ † φW a µν W a µν . Our notation for field operators will follow the notation in CKW. There is one figure for each operator consisting of up to six panels. The captions for each panel are the same for each figure and described below. In addition, we make operator-specific comments in some individual figure captions.

A. Panel A:
We show the non-relativistic DM annihilation cross sections σv as a function of DM mass for each of the annihilation modes allowed for the particular EFT operator, assuming that a) the DM annihilates only through that operator and b) the DM particle is a thermal relic with Ω DM h 2 = 0.12.
If all annihilation modes are pure s−wave, the total cross section for annihilation attains a value of around σv ≈ 3.6×10 −26 cm 3 s −1 for M = 100 GeV [showing logarithmic variation with M ; cf. Eq. (1)]. However, p−wave components to some annihilation channels can cause the the present-day total annihilation cross section to be suppressed due to their non-negligible impact in the early universe when the relic density is set; these contributions naturally drive Λ larger and hence the present-day annihilation cross section smaller. This effect is particularly pronounced in operators where the s−wave annihilation is essentially independent of M but the p−wave annihilation contribution grows with M , e.g., the operator XVII-1, Ignoring kinematic thresholds, for all final states involving fermions except those in Table IX of CKW, the fermion contributions separate into four distinct types: leptons with T 3 = ±1/2 and quarks with T 3 = ±1/2 (cf., the results for all operators XVIII-XX).

B. Panel B:
The per-annihilation total differential spectrum dN γ /dE γ = f BR f ·dN f γ /dE γ of prompt gamma rays as a function of E γ is plotted for different values of DM mass, indicated by the color scale. The mass dependence of the spectra depends on whether the dominant final states are fermions or if they are gauge bosons and Higgs. In the latter case, the spectrum becomes harder with increasing DM mass because the final states become more boosted. Further details on these spectra can be found in Sec. II D. No inverse Compton scattering (ICS) component is included. An estimate of the error introduced by neglect of ICS is given in Appendix A.

C. Panel C:
If applicable, this panel shows the line search limits for the DM mass region 5 GeV to approximately 300 GeV from the Fermi-LAT [16] line search using observed photon fluxes from the Galactic center. We show these limits in terms of the total annihilation rate to monochromatic photons (i.e., the annihilation rate to each final state weighted by the number of monochromatic photons in that final state). The solid black line shows the specific σv required for a thermal relic.
The solid colored lines show 95% CL UL for various halo profiles and ROIs, assuming the central values for halo normalization ρ from Ref. [31] (see Table I). The like-colored bands show the variation in the UL as ρ is varied through 2σ limits. For operators with multiple possible lines (i.e., γγ and γZ annihilation modes) we do not set limits in the region of DM masses 80 GeV to 160 GeV where we expect the Fermi-LAT line shape would become a very broad or double-peaked structure (see Sec. III A).
If present, the colored squares with their error bars represent (M, σv) values where the operator in question could additionally supply the line reported in Weniger [20], assuming the central values for ρ for each halo profile. As discussed in Sec. III F, for operators where there are multiple photon lines we make the assumption that the line at 130 GeV is due to the dominant photon line mode. If applicable, this panel shows the line search limits for the WIMP mass region from 500 GeV to 20 TeV from the H.E.S.S [17] line search discussed in Sec. III B. The ROI is a circle of radius 1 • centered on the GC with Galactic latitudes |b| < 0.3 • excluded. Note the change from Panel C to a logarithmic scale in M .
The black and solid colored lines and similarly colored bands are all as described in Panel C. Isothermal limits are very weak for this line search owing to the cored nature of the profile and very small ROI near the GC.

E. Panel E:
This panel shows the inclusive (continuum and line) spectrum limits for DM mass from 5 GeV to 1 TeV based on the "model-independent" limits given in Hooper et al. [18] and described in Sec. III C. These are given as the solid colored lines labelled 'GC ...'; the colored bands show the range of halo profile normalizations, as in Panel C. The solid black line is the result for the total σv for the operator.
The sharp features in the GC limits that can be seen in some cases are due to the γ lines changing through the bins used in [18]. For operators with lines only from γZ or γh, we also see sharp features in the limits near threshold as the line energy migrates rapidly across the bins; had we convolved the lines with the Fermi-LAT energy dispersion, these cross-over regions would be smoothed.
Also shown are Fermi-LAT [19] stacked dwarf galaxy limits, discussed in Sec. III D. These are given as a variety of grey lines and are a conservative estimate of the 95% CL UL limits; because the published limits were presented in terms of specific final states, there is not enough information to fully derive the limits for the continuum spectra here.

F. Panel F:
If applicable, the limits presented here are on the ratio of selected continuum final state annihilation cross sections to the line final states γγ and γZ as described in Cohen et al. [21] and summarized in Sec. III E; these limits are only applicable if the operator has no final state γh, and has dominant annihilation branching ratios to at least one of the final states W W/ZZ or bb.
The solid black line is the ratio R th ≡ (σ W W v + σ ZZ v + σ bb v) / (2σγγv + σ γZ v) for this operator; the dashed colored lines are individual annihilation mode contributions to this summed result. The various solid colored lines give either the supersaturation 95% CL UL on R for fixed M , or the boundary of the '2σ' confidence region 12 in the R − M plane for the shape constraints from Cohen et al. [21] assuming 100% annihilation to the indicated final state. Further discussion of the applicability of these limits is given in the text in Sec. III E, but it should be borne in mind that operators which are only marginally excluded or allowed per the results given here merit further investigation to determine whether or not they are actually excluded. In particular, the continuum radiation from the Z in γZ is not factored in. Furthermore, the ratio of γγ to γZ is not fixed at the correct ratio for each operator here considered, but is rather marginalized over in the analysis of Ref. [21].

V. DISCUSSION AND CONCLUSIONS
Broadly speaking, we find that the 34 operators with s-wave annihilation channels can be classified into thirteen qualitatively distinct classes on the basis of the limits which apply to each operator; we summarize these in Table III. • Operators coupling to hypercharge gauge bosons; photon lines present (categories 1-3) For operators in the first 3 categories in Table III, annihilation to final states giving lines has a large branching fraction as soon as the dark matter mass is larger than lowest kinematical threshold for such channels. This is a consequence of the fact that the B field strength tensor has a dominant photon component so that couplings to the Z are parametrically suppressed with respect to couplings to the γ by a factor of tan 2 θ W ∼ 0.3 (this suppression may be partially invalidated by other large factors arising in the evaluation of the matrix element, but nevertheless serves as a useful rule-of-thumb). As a result, the operators in categories 1 and 2 are excluded at 95% confidence by the Fermi-LAT [16] and H.E.S.S. [17] line limits for most masses from the lowest kinematic γX and ZX where X = Z (for "−h.c." operators) OR X = h (for "+h.c.") are present above thresholds, with the latter in each case typically suppressed by (sin θW / cos θW ) 2 . Line limits are constraining up to a TeV.
γX, ZX and W + W − where X = Z (" − h.c.") OR h (" + h.c.") are present above thresholds. Line limits are severely constraining for masses less than around 100 GeV; once the W + W − mode enters, the γX annihilation mode is heavily suppressed. Inclusive limits are weakly constraining for masses up to a few hundred GeV. Where applicable, the ratio limits are either marginal or quite constraining on 125-150 GeV.
Similar to the category just above, but p-wave components to the cross section cause suppression of the s-wave cross section and weaken line constraints above a few hundred GeV. threshold for a line final state up to a few TeV, except a) if the very conservative isothermal profile choice is made for the H.E.S.S. limits and b) the dark-matter mass is in the region 300 GeV and 500 GeV, which is not addressed by any of the line limits. Depending on the profile choice, the required thermal cross section for these operators may also be in tension with the experimental upper limits from a few to 20 TeV.
For operators in category 3, the total s-wave cross section drops away from the canonical thermal value of approximately 3 × 10 −26 cm 3 s −1 as M −2 for M above a few hundred GeV. This phenomenon occurs generically for operators with fermion DM coupling via an axial-vector J DM current, and is due to the presence of p-wave components to the annihilation cross section which are enhanced over the s-wave by a factor of s/m 2 V ∼ M 2 . This enhancement can be traced to a coupling in the p-wave (but not the s-wave) to the longitudinal component µ L ∼ k µ /m Z of the massive vector. Since both s-and p-wave components of σv can be important at freeze-out [cf. Eq. (1)], having a/b ∼ M −2 naturally suppresses the present-day annihilation cross section for large M by driving up the value of Λ required to saturate the present-day average DM density. As a result, a significantly smaller range of M is definitively excluded by, or put in tension with, the experimental limits compared to categories 1 and 2, especially for TeV-range DM. The constraints from the inclusive limits from Ref. [18] are less stringent: in most cases, the limits are in tension with the thermal relic cross section for M below one or two hundred GeV only for some choices of halo profile, except for operators with γγ modes where, below ca. 100 GeV, the inclusive limits exclude such operators even for conservative profile choices.
• Operators coupling to SU (2) L gauge bosons and Higgs; photon lines present (categories [6][7][8] For operators coupling to the SU (2) L gauge bosons and Higgs (entries 6-8 in Table III), the W + W − final state is always dominant for M m W since in this case couplings to the Z and γ are parametrically suppressed by cos 2 θ W and sin 2 θ W , respectively. Consequently, for operators in categories 6 and 7, the line limits are only severely constraining between the lowest kinematic threshold for any line mode and M m W ; for M greater than about 100 GeV and less than a few TeV the experimental limits are typically either only in tension with the required thermal cross section for more aggressive choices of halo profile, or are completely unconstraining. Operators in category 8 again contain fermion DM with coupling via an axial-vector current and thus suffer the same falling σv at large M as the operators in category 3, and as a result are even less constrained at large M than the operators in categories 6 and 7. The inclusive limits from Ref. [18] are somewhat stronger for all these operators than for the operators in categories 1-3, and may be in tension with the required thermal cross section up to a few hundred GeV depending on the profile choice. Nevertheless, taken together with the line limits, these operators are only weakly constrained for heavy DM.
In addition, for most operators in categories 6-8, the shape-based ratio limits [21] are either very constraining (operators XI-3 and XVII-3), or are marginal in the sense that the current limit is very close to the predicted value from the operators. We however caution the reader that the limits in Ref. [21] were obtained by assuming the line signal would account for the Fermi signal for some choice of DM mass, and allowing the relative strengths of the lines to float to the optimal value. The first assumption may or may not be applicable to our analysis depending on the operator considered; but the second is never applicable to our analysis since the relative strengths of all annihilation channels are fixed for any given operator. Therefore, marginal cases may or may not actually be constrained by these limits and a more careful analysis is needed here.
• Operators with annihilation modes to fermions; photon lines present (categories 4,5,9,10) The analysis in Ref. [18] as applied to our operators indicates that if the only available s-wave annihilation mode for M < 10 GeV dark matter is to Standard Model fermions, there is a strong constraint which seems to rule out the annihilation cross section necessary for obtaining the thermal relic abundance for M < 10 GeV.
The operators listed in entries 4, 5, 9 and 10 of Table III all have annihilation modes to both fermions and line final states (amongst others). Of these, the operators coupling to the field-strength tensor (rather than its dual) have s-wave annihilation to ff final states and so can be excluded for M below 10 GeV; they may also be in tension with the inclusive limits, depending on the profile choice, from 10 GeV to a few hundred GeV. Also, for these operators, the branching ratio to ff annihilation modes falls once the on-shell diboson modes become kinematically allowed. This is because the annihilation to fermions goes through s-channel gauge boson exchange and so requires an additional Higgs vev insertion and coupling of fermions to γ/Z compared to the Zh and γh modes (see the next paragraph for a comment on the W + W − mode). On dimensional grounds this results in a suppression of the ff modes relative to the Zh and γh modes by a factor ∼ m 2 Z,W /s (where s ∼ 4M 2 ). For the operators coupling instead to the dual field strength tensor, the fermion cross section is p-wave so the inclusive limits are less constraining: they may be in tension with the required cross section from the γh threshold to a few hundred GeV, but only for fairly aggressive profile choices.
For the operators in categories 4 and 5 with couplings to hypercharge and Higgs, the photon line becomes strong above the kinematic threshold for the γh annihilation mode leading to strongly excluding line constraints, whereas for the operators in categories 8 and 9 with couplings to the SU (2) L gauge bosons and Higgs, the diboson annihilation modes are dominated by W + W − , which leads to line limits for these operators in tension with the required thermal cross section only if an aggressive profile choice is made. We note that this occurs notwithstanding the fact that for the W + W − mode (as for the ff mode), both Higgs fields must be replaced by their vevs; the reason that this does not suppress this final state is that there is a cancellation of the suppression factor ∼ m 2 W /s against an enhancement ∼ s/m 2 W from the longitudinal mode of the W ± (the by-now-familiar suppression, relative to W ± , of couplings to Z and γ by weak-mixing-angle factors accounts for the relative strength of the W + W − and Zh or γh modes).
• Operators without photon lines (categories [11][12][13] Finally, entries 11-13 in Table III are operators which couple only to SM final states without photon lines. 13 The two operators in entry 11 have strong annihilation to ff either for all M (operator XX-1) 14 or at low M (operator IX-2), 15 but in either case this means that they are strongly constrained for M 10 GeV even for conservative profile choices. Owing to the very flat (with changing M ) nature of the inclusive upper limits on σv, the astrophysical uncertainty gives a rather broad range of M for which the operators may be in tension with the experimental limits depending on the profile choice: 10 GeV to a few hundred GeV. The operators in entries 12 and 13 have only Zh final states: the former may be in tension with the astrophysical limits for aggressive profile choices for M between 100 and 400 GeV, but they are unconstrained for more modest profile choices; the latter are unconstrained regardless of the profile chosen due to s-wave cross section suppression by the same enhanced-p-wave effect which was discussed at length above.
In addition to these indirect limits, we emphasize the point made in CKW that the 'tensor' fermion DM currents χγ µν χ (e.g., operators XVIII-n, XIX-n and XX-n in categories 4-5 and 9-11) give rise to magnetic or electric dipole couplings (if necessary, by replacing both Higgs fields with their vevs) which causes them to be very strongly constrained by direct detection experiments [42,43]. Given our assumptions about the DM being a cold thermal relic WIMP saturating the measured average density, the operators giving electric-dipole couplings (coupling to the dual field tensor) are excluded at at least 90% confidence by both CDMS and XENON across the entire mass range from M = 10 GeV to at least 20 TeV (probably somewhat larger), while the operators giving magnetic-dipole couplings (coupling to the field tensor, not its dual) are excluded at at least 90% confidence from M = 10 GeV to M > 20 TeV, M 1 TeV, or M 400 GeV for operators XX-1, XVIII-1, or XIX-1, respectively (with the exception of a very narrow mass window around m Z /2 for the latter two of these).
We have presented all of our results under the assumption that fermionic DM is Dirac, however the compensating factors of 2 in Eqs. (1) and (6a) imply that there is a uniform shifting of all σv values down by a factor of 2 if Majorana fermion DM was assumed instead; therefore, the same thermal cross section relative to all of our 95% CL UL limits would obtain, at least for all the operators which do not vanish identically for Majorana fermion DM.

γ-ray Line Signals
As discussed in Sec. III F and Sec. III G, there is some evidence, first reported by Weniger [20], for a photon line in the Fermi-LAT data near 130 GeV (updated to approximately 133 GeV in Ref. [16]) with a flux matching DM annihilation into monochromatic photons with cross section σv ≈ 10 −27 cm 3 /s. Although the most recent official Fermi-LAT collaboration analysis [16] finds less significant evidence for this putative line-like feature, it cannot fully explain it as a systematic effect, and cautions that more work is needed to understand it. If we however take this signal seriously, we find that there are a number of operators (indicated in Table III  The signal strength for the photon line is a factor of a few (up to an order of magnitude, depending on the halo profile) below the thermal relic cross section. Operators with annihilation into SU (2) L gauge bosons (categories 6-8) naturally give a suppression of this size simply from the mixing angle sin 2 θ W . Numerical factors are also important, and those operators with smaller branching ratios into lines do not produce enough monochromatic photons if the cross section is fixed by the matching onto the correct relic density. To put it another way, these operators may be in 13 Annihilation to monochromatic photons is also possible at subleading order, for example through loops or see Ref. [23]; we do not include these subdominant modes. 14 Here we include annihilation modes purely through s-channel gauge boson exchange which explains why the suppression of the ff modes does not occur for large M as it does for operators XVIII-1, XIX-1 and IX-2. 15 The explanation here for the suppression of the ff branching ratio for M ∈ [m W , mt] or for M mt is that this annihilation mode is via s-channel Higgs exchange which leads to a suppression from the hff coupling ∼ m 2 f /s. some tension with continuum limits or the constraint on the continuum-to-line ratio since the BR to such continuum modes is then larger. Operators XV-3 and XV-4 are most promising as explanation of the line, but for completeness we also note that almost any of the other operators in categories 6-8 in Table III could also work, including VI-3, VI-4, VIII-3, VIII-4, XI-4, XIV-3, XIV-4, and XVII-4. The operators XI-3 and XVII-3 have smaller branching ratios to lines, less than 10%, and so do not fit the line as well for the reasons discussed above. However, given the systematic uncertainties, it is not yet conclusive that these operators cannot explain the Fermi gamma-ray line. (Finally, although the operators XVIII-1 and XIX-2 also have lines of roughly the correct strength, they are ruled out by direct detection experiments, as discussed above.)

Other Indirect Constraints and Outlook
Photon fluxes do not of course supply the only prospects for indirect detection of dark matter. At present, the IceCube experiment [44,45] reports limits on σv from neutrino fluxes which are at best on the order of 10 −22 − 10 −23 cm 3 s −1 , which are thus much too weak to be constraining. The absence of sharp spectral features in the AMS-02 positron fraction data provides a constraint [46] on annihilation modes giving a sharp edge in the positron fraction: we expect that these limits could be stronger than the GC inclusive limits for M less than a few tens of GeV for operators with sizeable BR for annihilation to e + e − , and to a lesser extent µ + µ − , final states (operators XVIII-1, XIX-1, IX-2 and XX-1; the latter may be more constrained since ff dominates for all masses). A slightly different analysis [47] constraining only the total number of positrons from various DM annihilation modes using a combination of Fermi-LAT and AMS-02 data confirms the approximate equivalence of the limits from positron measurements and the GC inclusive photon flux limits for the µ + µ − annihilation mode, but indicates that constraints arising from positron measurements for other channels (bb, W + W − , ZZ) would be between one and two orders of magnitude weaker than the GC inclusive photon flux limits. Additionally, constraints from antiproton fluxes may be useful to consider. Ref. [22] found that the PAMELA antiproton limits are typically weaker than the inclusive GC photon flux limits for W + W − or bb annihilation modes and do not typically have sufficient reach to exclude the thermal WIMP scenario for at least some of the operators we have considered; however, Ref. [48] finds that antiproton bounds are roughly comparable with the GC inclusive photon limits from Ref. [18] for pure bb/tt or W + W − /ZZ final states for M 100 GeV (the antiproton limits are much weaker for smaller M or for leptonic annihilation modes), but which have fairly large uncertainties arising from cosmic ray propagation models. One would have to fully reconstruct the analysis of Ref. [48] to take advantage of these comparable limits though, since independently comparing the individual mode cross sections to their respective limits would again be subject to the cost of a factor of BR f , which is usually fairly small for bb or tt modes for our operators.
Looking forward, Ref. [49] indicates that it would be appropriate to consider a factor of 3-5 improvement in the present best limits from the next generation of indirect detection experiments looking at photon fluxes (GAMMA-400, CTA, H.E.S.S.-II). This magnitude of improvement would be able to constrain a fair amount more of the parameter space for the mass of the thermal relic WIMP scenario even factoring in the astrophysical uncertainties, but would still not be able to completely rule out all DM masses for every EFT operator with s-wave annihilation channels which we have considered.                                    Λ −4 (φ † ∂ µ φ + h.c.)(W a λµ H † t a D λ H + h.c.) -:-Operator X -2                                                    Figure captions are provided in Sec. IV A through Sec. IV F. This operator may be compatible with present experimental limits and account well for the 130 GeV photon line for some profile choices; however, this operator gives rise to a magnetic dipole coupling of the DM and as such is excluded at 90% confidence by the CDMS and XENON direct detection experiments [42,43] for M from 10 GeV to ca. 1 TeV.     Figure captions are provided in Sec. IV A through Sec. IV F. This operator gives rise to a magnetic dipole coupling of the DM and as such is excluded at 90% confidence by the CDMS and XENON direct detection experiments [42,43] for M from 10GeV to ca. 400GeV.   Figure captions are provided in Sec. IV A through Sec. IV F. This operator may be compatible with present experimental limits and capable of accounting for the 130 GeV photon line for some profile choices; however, this operator gives rise to an electric dipole coupling of the DM and as such is excluded at 90% confidence by the CDMS and XENON direct detection experiments [42,43] for the entire range of M shown here.  of the ICS to prompt signal when considered over O(10 • ) regions. We therefore take the approximation of the limits in Ref. [18] being applicable to the combined ICS and prompt photon signals.
Since we are not however reproducing the template subtraction analysis, we must make a choice for the ROI over which to extract the effective ICS spectrum in order to make a direct comparison with the cross section limits in Ref. [18]. Motivated by the cuspy nature of the expected signal morphologies, and the fact that the photon signal used in setting the limits in that reference is mostly contained within a region of about 3 • around the GC, we chose to utilize a 5 • radius circular region centered on the GC in extracting the effective ICS spectra. Furthermore, at very high injection energies, the high-energy ICS signal is more constrained to the galactic center, since electrons must be at high energy to create high energy photons via single Klein-Nishina-regime Compton scatters, and this only happens near the point of production as the electrons lose energy as they propagate out of the GC region.
We integrate the ICS signal over this ROI and normalize by the line of sight integral for that region:  95% CL UL (so f < 0 means a more stringent limit with the effect properly accounted for, while f > 0 means the limit as currently set is too stringent). Values estimated by us are marked with a ; values taken from the relevant literature source are not marked. For comparison, the typical variation between the most conservative normalization for the most conservative profile choice and the most aggressive normalization for the cuspiest profile is about an order of magnitude, so relative to the mid-point the astrophysical uncertainties can cause the limit to vary by about +100% / −75%.

Result(s) affected
Description of effect Variation in result f To check our work, we compare with the analysis of a single line presented in Ref. [20]. We compute the value 21 of σ γγ v = (8πM 2 /2J)N γγ /E ROI assuming θ γz/γγ = 0, N γγ = 33 and M = 130 GeV, along with the Einasto profile J factor at the bottom of Table V which has been computed using the indicated parameter values taken from Ref. [20]. We obtain σ γγ v = 3.2 × 10 −27 cm 3 s −1 , which is a factor of 3 larger than the central value (σ γγ v = 1.1 × 10 −27 cm 3 s −1 ) reported in Ref. [20]. There is thus some tension between these two analyses.