{'fraction_non_alphanumeric': 0.030864197530864196, 'fraction_numerical': 0.014980427582053599, 'mean_word_length': 4.342941413058051, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improved Constraints on Dark Matter Annihilation to a Line using Fermi-LAT observations of Galaxy Clusters

Galaxy clusters are dominated by dark matter, and may have a larger proportion of surviving substructure than, e.g, field galaxies. Due to the presence of galaxy clusters in relative proximity and their high dark matter content, they are promising targets for the indirect detection of dark matter via Gamma-rays. Indeed, dedicated studies of sets of up to 100 clusters have been made previously, so far with no clear indication of a dark matter signal. Here we report on Gamma-ray observations of some 26,000 galaxy clusters based on Pass-7 Fermi Large Area Telescope (LAT) data, with clusters selected from the Tully 2MASS Groups catalog. None of these clusters is significantly detected in Gamma-rays, and we present Gamma-ray flux upper limits between 20 GeV and 500 GeV. We estimate the dark matter content of each of the clusters in these catalogs, and constrain the dark matter annihilation cross section, by analyzing Fermi-LAT data from the directions of the clusters. We set some of the tightest cluster-based constraints to date on the annihilation of dark matter particles with masses between 20 GeV and 500 GeV for annihilation to a gamma-ray line. Our cluster based constraints are not yet as strong as bounds placed using the Galactic Center, although an uncertainty still exists regarding the"boost factor"from cluster substructure, where we have chosen a rather conservative value. Our analysis, given this choice of possible boost, is not yet sensitive enough to fully rule out typical realistic DM candidates, especially if the gamma-ray line is not a dominant annihilation mode.


I. INTRODUCTION
Weakly Interacting Massive Particles (WIMPs) (for reviews, see Refs. [1][2][3]) are thought to be among the best motivated dark matter (DM) candidates. Some of these candidates annihilate among themselves in the early universe with a strength to naturally provide the correct relic density today to explain the dark matter abundance in the universe. It is natural to assume then that these DM candidates go through similar annihilation processes in the present universe, wherever the DM density is sufficiently high. This is the basis for DM indirect detection experiments, which search for the annihilation products of such dark matter particle candidates. Possible interesting annihilation model final end products include, but are not limited to: electrons and/or positrons, antiprotons, photons, and neutrinos.
Other than the Galactic Center, the lack of excess signal over background is used to place γ-ray flux upper limits for energies between 500 MeV and 500 GeV. These bounds are then used to constrain the dark matter annihilation for a broad range of particle masses and annihilation channels. For current bounds on DM annihilation, see the Fermi-LAT combined analysis of dwarf galaxies [54], VERITAS [55], and MAGIC [60] observations of Segue 1, as well as results of the H.E.S.S. collaboration [61] on Sagittarius and other dwarf galaxies. In addition, Ref. [56] showed that Fermi/LAT observations of Dwarf Galaxies highly constrain a dark matter interpretation of excess positrons seen in AMS-02 [57], HEAT [58], and PAMELA data [59].
In this paper, we focus on searches for γ-ray signatures of DM annihilation from galaxy clusters in Fermi-LAT data. These include the most massive virialized DM structures in the universe and could produce substantial DM annihilation signals. We use a dark matter annihilation model within galaxy clusters from which luminosity is derived from two components. The "smooth" component can be described by a radial model within each cluster. The "clumpy" component is described by a large number of subhaloes predicted by numerical simulations to exist within the clusters. The annihilation signal is likely enhanced by the "clumpy" component from the substructure, and this is expected to dominate over the smooth component by a large enhancement, or "boost" factor. Previous searches finding no evidence for gamma-rays from clusters have been used to set bounds on DM annihilation [62][63][64][65]. A novel feature in the present analysis is that, whereas previous work examined at most ∼ 100 clusters in total, we enlarge the sample to encompass tens of thousands of galaxy clusters. We also investigate directly these stacked data for structure in the energy distribution from the clusters (which all are at low redshift enough to have negligible energy shift compared to the ∼10-20% instrument resolution). We note that our strategy is different from that of recent previous studies, where attention is instead given to angular correlations between galaxy gamma-ray data and data of other wavelengths [66,67]. It is more similar in spirit to a recent analysis of the Fermi-LAT collaboration where, however, only 16 nearby clusters were studied [68]. In fact, one of the first analyses of this kind [69] claimed observation of a gamma-ray line signal around 130 GeV from galaxy clusters that would have been consistent with the indications from [70,71]. These interesting indications have, however, not been verified with further Fermi-LAT data [72].
The motivation for our analysis is the large number of galaxy groups identified by Tully in the 2MASS Redshift Survey (2MRS) data [73]. We search in the direction of these tens of thousands of galaxy groups (hereafter denoted clusters) for an excess of γ-rays over expected backgrounds. In this first paper, we look for a γ-ray line [13,14] where, due to the low velocity of the annihilating particles, the energy of the photons is equal to that of the WIMP mass, In a forthcoming paper, we will generalize to searching for gamma-rays from DM annihilation through a variety of other channels in addition to annihilation directly to a line. For the IR-selected galaxies in the Two Micron All Sky Survey (2MASS), redshifts are provided by the 2MASS Redshift Survey (2MRS) [74]. The 2MRS galaxies have been selected to be at Galactic latitude |b| ≥ 5 • (≥ 8 • towards the Bulge). Previous work studying DM annihilation with 2MASS data includes [75]. Using a combination of 2-D spatial position from 2MASS together with the redshifts, Tully was able to identify galaxies grouped together in clusters [73]. Given their redshifts, one knows the distances to the clusters. Errors on the distance will be related to redshift space distortions (bulk movement of galaxies) and the uncertainty on the Hubble constant. These clusters all have redshifts consistent with being nearby, within 280 Mpc (most within 50-100 Mpc). The IR light from galaxies is used to determine their stellar masses, from which one can obtain the total mass using abundance matching. This mass estimation is quite accurate because of the lack of extinction in IR. For the largest masses, velocity dispersion can be used to estimate the total mass of the cluster, and for smaller groups more detailed kinematics of nearby groups can be used to find empirical estimates [73]. The estimates should be accurate to ∼ 20%. The clusters in Tully's catalogue have masses in the range 10 12 − 10 15 M . With the halo masses and distances provided in [73], we will thus be able to calculate the predicted γ-ray flux from dark matter (for a given WIMP mass and annihilation cross-section).
Our general approach is as follows. For each cluster, we define a region of interest containing 95% of the expected DM annihilation luminosity. We slide a bin of energy range ∼ twice the energy resolution of Fermi-LAT across our full spectrum of interest, centered at photon energies from 20 GeV up to 500 GeV in steps of 2 GeV. We then add up all the observed photons for all the clusters for each energy window, and compare our observed photon count to the expected count. We search for a bump in the observed photon count above an expected power law background, i.e. we search for a line (or internal bremsstrahlung [70,76]) signal at an energy equal to the WIMP mass due to DM annihilation. Since no excess above background is found, we use the null signal to place bounds on the DM annihilation cross section as a function of WIMP mass.

II. DARK MATTER ANNIHILATION γ-RAY FLUX
The differential γ-ray flux dφ k dEdΩ from DM annihilation in cluster k can be written as the product of two components, a factor dΦpp dE that encodes all the particle physics and the so-called J-factor [13] that contains the astrophysics, Here, σv is the DM annihilation cross section, m χ is the DM mass, and the δ function indicates that we are interested in a line signal. The J-factor (the term in the second set of parentheses) integrates the square of the DM density ρ χ along the line of sight. For each cluster we must now obtain an estimate of the total J-factor. The J-factor has two contributions, the smooth halo component and the substructure (clumpy) component.

The Smooth Component
From simulations [65,80], for the smooth component of our DM haloes we take the Navarro, Frenk, and White (NFW) density profile where the characteristic density ρ s and radius r s are related to halo concentration and virial radius through the relations, ρ s = 200 3 and r s = r 200 /c. Here ρ c is the critical density of the universe, r 200 is the cluster virial radius within which the average density is 200ρ c , and the concentration parameter c is given by [82]. Here M 200 is the virial mass given by Tully.
The J-factor integrated over a large enough solid angle can be used to define the total flux, J int = ∆Ω JdΩ (this value would be appropriate for a point source approximation). For the smooth component of a cluster, Han et al. [65] find that (see also [16]) Here D A is the angular diameter distance of the cluster.
J N F W is a quantity we will use for cluster selection and normalization for the boost factor described shortly.

The Clumpy Component
We must now account for enhancement due to substructure inside the cluster, which is expected to account for the majority of the dark matter signal from galaxy clusters [17][18][19]65]. While the total sub-halo mass constitutes only 10 to 20 percent of the total halo mass, the DM density inside the subhaloes is considerably higher density than within the smooth component of the cluster. Indeed, in comparing with data, it is the substructure J sub that is the dominant component in the last term of Eqn. (2).
We define a boost factor as the ratio of the total integrated luminosity due to the substructure (subscript sub) to the total integrated luminosity due to the smooth NFW component. The notation is as above, with the total integrated flux from the substructure (over a large enough solid angle) as In the last few years much progress has been made in determining the boost factor of dark matter annihilation in clusters (see [81] and references therein). The improvements have been a result of improvements in estimating the number of substructures in a halo, and the annihilation rate in a give sub-halo(which depends strongly on the sub-halo concentration parameter) [18,19,65]. The best-guess boost factor per halo is between about 20 and 50 depending on the mass of the halo in our mass range, with boost increasing with cluster mass. Thus in our work we use the conservative boost factor of [81,84] ranging from B=20 for 10 12 M clusters to B=55 for 10 15 M clusters (as opposed to the less conservative one of Bi et al. [85] which is roughly a factor of 20 or more larger).
As substructures are to a large extent destroyed by tidal effects near the central parts of the clusters, the DM photon annihilation profile will be substantially more extended than the NFW mass profile [18,19,65]. We are therefore using the lack of observed line signals in a relatively large region around each cluster to place bounds on the DM contribution. We consider a region of interest which will provide at least 95 percent of photons from cluster which as mentioned typically corresponds to about a degree on the sky.
Simulations of [83] show that the surface brightness profile of sub-halo emission can be fitted within r 200 by Following [65], we take for the subhalo emission surface brightness beyond the virial radius an exponential decay, The Total Differential Flux and Cutoff Radius The total annihilation profile, which is the sum of the contributions from a smooth NFW profile and the subhalo emission, is completely dominated by subhalo emission except in the very center of the cluster where the smooth NFW cusp is important.
Rather than accumulating all the photons emitted by a cluster, we perform a signal cut at a radius from the center of the cluster that contains 95% of the DM signal (from the substructure). Thus we take the expected differential flux of photons from a cluster k due to DM annihilation to a line to be where is the integrated J-factor that includes 95% of the photons from the substructure of the cluster. Given J 0.95 from theoretical estimates of DM annihilation, we now look for a conical section from the point of view of the Earth, defined by and angle θ max , that satisfies This angle then defines a radius that contains 95% of the DM annihilation signal from substructure (We are not including instrumental error in r max . We will return to the effect of the point spread function when defining the region of interest). The angular diameter distance to the cluster D A is obtained from data as described in the next section. We will cutoff any signal beyond this radius r max , which is typically roughly twice the viral radius of the cluster. In our analysis, we ignore the smooth component of the halo, thereby underpredicting the DM signal and obtaining slightly more conservative bounds. 1

III. DATA SELECTION
We perform what is denoted as a stacked cluster approach in our data selection and analysis. Previous work done using clusters to constrain dark matter annihilation into a line have taken two different approaches: a detailed cluster approach or a stacked cluster approach. For a given set of clusters, the detailed approach has stronger bounds than the stacking method. In the detailed case, each halo is carefully analyzed and finally a joint likelihood is performed on tens of objects [77]. These authors have typically used galaxy clusters which are bright in X-ray (also presumably also dark matter rich). In the stacking approach, the data is stacked before performing joint likelihood analysis. One such stacking approach is [78], which looks at a cross correlation between photons from the Fermi-LAT with galaxy surveys such as 2MASS and SDSS.

A. Selecting Galaxy Clusters
We use Tully's 2MRS catalog of clusters or groups of galaxies, which is widely regarded as one of the best and most comprehensive catalogs of galaxy clusters derived from 2MASS data. 2MRS includes redshift information. Tully's catalog of 2MRS data is optimal for our work for two reasons. First, the large number of clusters guarantees that we will find promising targets from the point of view of large dark matter annihilation signal and no coincidence with background point sources as described further below. Second, in a proper analysis we would like to have detailed and unbiased statistical description of the halo properties which will play a part in the dark matter signal, namely the halo's mass and distance. Some authors have previously exploited the 2MASS data (for example [79]) to construct galaxy cluster catalogs; however, Tully obtains improved estimates of the mass and distance to the clusters (see [73] for more details on the various statistical tests used in his work). We also note that 2MASS data have advantages over other surveys previously used to do studies of DM annihilation from clusters -clusters in X-ray catalogs, SDSS data, and using the S-Z effect. 2MASS has the two advantages of providing a larger number of clusters and, since they are nearby (z < 0.01 ), they have higher luminosities. For example, the 10 brightest clusters in our study are nearly three times brighter than the X-ray clusters used in the work of Anderson etal [68]. Here brightest means the highest J N F W . Table III A shows the five brightest clusters used in our analysis. Indeed none of our 16 brightest clusters overlaps with the 16 clusters used in Anderson etal, as shown in Fig. 4 in comparison with Fig. 3. One reason for the lack of overlap is that 2MASS contains so many clusters that we simply throw away those that overlap with any of the 150 brightest objects in the Fermi/LAT point source catalog. Whereas two of the brightest objects in the sky from the point of view of DM annihilation are Virgo and Fornax, we do not include them in our studies due to the point source contamination. Anderson etal dealt with the point sources by doing careful modeling of the clusters and cutting out the point sources using spatial information. We, on the other hand, throw them out completely in favor of using some of the many other clusters in the 2MASS data. In principle it would be interesting to add them back in to obtain tighter bounds on DM annihilation from the combination of the data.
Henceforth our studies are performed for three different cuts on halos: all the halos in the Tully 2MASS catalog, the 2500 brightest halos, and the 600 brightest halos. Again, we use the word "brightest" to mean those clusters with the highest value of J N F W . Table II illustrates the total value of J N F W added up for all the clusters in each of the three cuts.      Fig.3, we consider the case where we also include clusters which overlap with one of the 150 point sources. With the same cut off as in Fig. 3, we now have 100 more galaxy clusters. Note that the brightest clusters are Virgo and Fornax in that order. In our analysis, we do not include Virgo and Fornax since they overlap with one of the bright 150 point sources.

B. Selecting γ − rays from Fermi-LAT
The LAT instrument aboard the Fermi satellite is a pair-conversion telescope measuring γ-rays in the energy range from 20 MeV to > 500 GeV. For a more detailed description, see [86,87]. We analyze five years of public Fermi-LAT Pass 7 reprocessed data taken between 4 Aug 2008 to 8 Mar 2012 in the energy range between 20 and 500 GeV.
We apply the celestial zenith-angle cut θ < 100 • in order to avoid contamination with the earth albedo, as well as the recommended quality-filter cut DAT A QU AL == 1. We use the U LT RACLEAN events selection, i.e. the highest quality data, from P 7REP U LT RACLEAN V 15. We use both front-and back-converted events. The selection of As mentioned previously, we sub-select again to reject background by removing the 5% brightest point sources identified by Fermi-LAT which are coincident with the clusters we are studying. Consequently we only include clusters which do not overlap with any of brightest 150 Fermi-LAT point sources. By doing so, we exclude some of the brightest clusters (from the point of view of γ-ray production due to DM annihilation) on the sky such as Fornax and Virgo, see Fig. 3. Still after this sub-selection, we retain the vast majority of the clusters in the catalogs included; on the order of a few hundred clusters were removed out of 26,000 we work with. Similarly, when we apply the cut down to only the 600 brightest clusters, we have lost about 100 clusters. Compare Fig. 3 with Fig. 4.

Region of Interest
To compare the expected number of photons to observations, we must determine which photons in the Fermi-LAT data to include in our analysis of a given cluster. There are two ingredients in this choice. Previously we found (see Eqn. (18)) the maximum radius r max from the center of the cluster that satisfies the requirement of retaining 95% of the signal from the substructure. Second, we must account for the instrumental smearing of the directionality of the signal due to the point spread function (PSF). We note that the Fermi-LAT PSF varies as a function of energy, and above energies of 40 GeV the PSF varies slowly.
To calculate the selection radius of photons for a given cluster, the most accurate approach would be to convolve the point-spread function (PSF) of Fermi-LAT with the DM profile of the cluster. We can, however, take a simpler approximation. After extracting the PSF from the Fermi Tools, for each value of the DM mass, we calculated the radius within which 95% of the photons lie due to the uncertainty of the detector r P SF,0.95 . Specifically, the 95 percent PSF containment angle is radial angle describing a circle around a point source such that 95 percent of the photons from that point source are observed in that angle accounting for instrument error. The PSF containment  angle typically ends up being on the order of a degree for energies above 100 GeV. Then for each cluster, we consider a Region of Interest (ROI) for incoming photons taken to be the sum Our region of interest corresponds to a radius within which we look in the Fermi-LAT data for photons from any given cluster in the Tully data.

C. Photon Number Count and Energy Spectrum
The photon number count is shown in Table III for a variety of different selection cuts for photon energies in the range 10 GeV < E γ < 1 TeV. Only photons from outside the Galactic Plane with Galactic latitude |b| > 10 • are included in all cases. The total photon count outside the Plane is shown, followed by photons associated with all clusters, photons associated with the 2500 brightest clusters, the 600 brightest clusters, and the 150 brightest point sources in the Fermi-LAT data. Finally, the photon number count from outside the Galactic Plane but unassociated with clusters is shown as well.   Table. The green curve illustrates the spectrum for all the Fermi-LAT observed photons. Black is for photons from outside the Galactic Plane( |b| > 10 • ) that are associated with galaxy clusters. Blue is for all photons outside the Galactic Plane (the sum of associated and unassociated). Red is for photons again outside the Galactic Plane but unassociated with galaxy clusters and unassociated with any of the 150 brightest point sources. We have normalized the curves such that the area under each curve will equal one. All of the curves have similar shapes except that the green curve appears to fall less steeply with energy, because it includes sources in the plane of the Galaxy. The large error bars in the unassociated spectrum (red curve) at high energy are due to the low number of photons in that case (as shown in Table III). Fig. 7 can be used to illustrate two points. First, within our energy windows, the spectra can roughly be fit by power laws. Thus it is reasonable to use a series of power laws (within each energy windows) to approximate the background. The unassociated photons in particular are a separate sample from the clusters we are studying and provide a reasonable background estimate. Our maximum likelihood approach to handle background will be described in the next section. Second, a line signal would show up as a bump in this figure and, at least by eye, there doesn't appear to be any. Our analysis will confirm that no line is in the data at any great significance, hence allowing us to place limits on DM properties.
We also illustrate how the spectrum will vary for our different cuts on the clusters and galaxies in Figure 8. We plot the three different case considered in the paper: the spectrum for all of the clusters, the spectrum from 2500 clusters, and for 600 clusters. The spectral shape appears to be roughly the same for the different cases, yet with fluctuations between the different data sets that will affect our resulting bounds.

IV. EXPECTED NUMBER OF PHOTONS FROM SIGNAL AND BACKGROUND IN THE FERMI-LAT
We follow the approach of Weniger [71], Section 2.5 to calculate the total number of expected events ν w in the energy window (w) between where we have summed over all of the clusters k and integrated over the solid angle dΩ k over each cluster.
Here we have defined: θ = impact angle of photon events with respect to the instrument axis, (dT /dθ) k = the observational time at this impact angle for a cluster k, A(E , θ, j) = the effective area, E = the true energy, j = f, b = denotes front-and back-converted events, E = the reconstructed energy D(E|E , θ, j) = the energy dispersion of the LAT. (16) We have also defined the following flux quantities: (I mχ ) k = differential flux from dark matter into a line (see Eq. 2 ) For a given dark matter mass, we picked a window which is sufficiently large to effectively model the background. The energy resolution of Fermi-LAT is roughly 10-20 percent over the energies we are considering (see Fig. 9 ). We take the energy window to be m χ /1.4 < E < 1.4m χ , but with a maximum value for the lower limit of the energy window E − max = 200GeV (due to a dearth of very high energy photons). With a 100 GeV line, our window would then go from 71 to 140 GeV. We have found that taking a larger or slight smaller window does not affect the overall constraint, as has been similarly shown in [71].
We average over the instrument angle and back-and-front converted events. Again we follow Weniger and define an effective exposure for a cluster k as and an effective energy dispersion for cluster k as For cluster k, the expected distribution of the conditional parameters θ and j for events with energy E is From the publicly available Fermi Science Tools 3 , we obtained the dispersion function D(E|E , θ, j), effective area A(E , θ, j) , and dT /dθ for each cluster. Now that we have obtained the effective energy dispersion function for each cluster in Eq.(19), we can obtain an "average effective energy dispersion function" D ef f (E|E ) ave over the entire sky. We first average over the clusters to find (dT /dθ) ave . With this new quantity, we then average over the instrument angle and back-and-front converted events to obtain the average effective dispersion function over the clusters: where and In Fig. 9, we show D ef f (E|E ) ave for the case (E = m χ = 100 GeV). Weniger argues that the average effective dispersion function is a good approximation at any location in the sky. Thus, for a given dark mass m χ , for all of the clusters we will use the same average dispersion function. We now calculate the expected number of photons in some energy window (w), Here the solid angle integral dΩ k is over the ROI for the given cluster. The effective exposure varies slowly across the sky; we note that we do not use the average effective exposure here, but instead the local one. Since the ROI of a halo is typically only of order a degree on the sky, we can treat the exposure as constant across a cluster and take it out of the integral over dΩ k . We take the differential photon number count due to the background component to be a power law (within an energy window ranging from E − w = m χ /1.4 to E + w = 1.4m χ ) of the form Different WIMP masses have different power laws for their background fits. Now for the DM signal, the differential photon number count within the same energy window would be The first term on the RHS accounts for the instrument and the second term accounts for the particle physics. The last term defines a new quantity (which we call the exposed integrated J-factor) that accounts for both the astrophysics and exposure for each cluster, Here we have used Eq. (10) to determine the integrated J-factor (J 0.95 ) k for cluster k. The notation (J X ) T is intended to guide the eye of the reader to the fact that this quantity has been obtained by summing the product of two quantities: the integrated J-factor for each cluster times its exposure. The subscript T refers to the fact that this a total quantity. In addition, we use J k X k = (J X ) k as a shorthand for the individual J-factor and exposure for a given galaxy k, which will be useful in describing the J-factor Likelihood. We now have our expected photon energy spectrum, The first term gives the total flux from dark matter annihilation from all of the clusters used in our analysis. The second term account for the observed background.

V. LIKELIHOOD ANALYSIS
We performed an un-binned analysis with our likelihood given as where N is the total number of observed photons in the energy window w between E − w and E + w , ν w is the number of expected photons in the window w, E j is the observed energy of a photon, dE j /dE is given by Eq. (28). The product is over all photons in w for a given DM mass m χ . The first term on the RHS (a Poisson distribution) gives the probability of observing N photons with an expectation of ν w photons. The second term accounts for the probability of observing any one of the N photons. The last term (L J X ) is the likelihood function for the exposed integrated J-factor given in Eq. (27), which we now derive below.

A. J-Factor Likelihood Function
We take a Gaussian probability distribution for the Exposed Integrated J-factor likelihood function where σ T incorporates the individual errors on the exposed integrated J-factor for each cluster and any other systematic errors. T refers to total as in total error for σ T . (J X ) is the free variable and (J X ) T is the expected integrated J-factor given in Eq. (27).
Formally the relative error for the exposed integrated J-factor for a cluster will go like where δ refers to the error for a given quantity. We have ignored the errors associated with the galaxy's exposure X k . We note that the J-factor for a given cluster goes like which is derivable from Eq. (5) by simply rewriting r s and ρ s in terms of the mass of the halo (M ). D is the distance to the halo. The second step follows from the Hubble law, where V is the velocity of the galaxy and H 0 is the Hubble constant. In deriving the above expression, we have neglected two effects: one-the dependence of the denominator of ρ s on M and two-the mass dependence of the boost factor on M . We have two arguments to justify our simplification. First, the boost factor is a very slowly varying function of M . Also the denominator is very weakly dependent upon mass (a very small fractional power). Secondly, error due to velocity and error due to the Hubble constant are larger. We first provide an expression of errors for galaxy k with We add the relative errors in mass δM/M and in velocity δV /V in quadrature since they are independent. We note that for most of the clusters (δM/M ) k ∼ 0.20 and (δV /V ) k ∼ 0.10, which we will justify below. Including the uncertainty in the Hubble constant, we find the total error, where we sum over all N of the clusters in our analysis. For smaller (larger) values of H 0 , we can see from Eq. (32) that the J-factors of the clusters increase (or decrease), so that they appear brighter (dimmer). Again, we add the error in H 0 in quadrature. The errors in velocity and halo mass are independent of the Hubble constant. In our code we numerically calculate Eq. (34). We find that the errors associated with H 0 dominate in σ T over those related to the properties of the cluster i.e. δV /V and δM/M . We now provide a heuristic argument that provides some intuition as to why the systematic error in H 0 leads to the largest contribution to σ T . As a starting point, we write the first term of Eq. (34) as N k σ 2 k = N J X 2 δJ /J 2 . Here we use the notation ... to indicate quantities averaged over all the clusters, e.g. the average integrated exposed J-factor is J X = (J X ) T /N . With appropriate substitutions, we can now write For N 1, the second term dominates. In our case, we have hundreds to tens of thousands of clusters (depending on the selection cuts), so that this criterion is easily satisfied. Errors on the Hubble constant are on the few percent level. With Eq. (33), we can estimate δJ/J which depends upon δV /V and δM/M . δM/M is on the order of twenty percent for the mass method of the Tully 2MASS Catalog. In the case of the velocity, δV /V is on the order of 10 percent 4

B. Putting Pieces Together
Now that we have explained the various elements which contribute to our overall likelihood function. We can now turn to constraining the properties of dark matter. Using a delta-log-likelihood approach we can use Fermi-LAT data to constrain σv for a given m χ by treating β, γ, and (J X ) as nuisance parameters, which we profile out. The delta-log-likelihood ∆L is given by where σv , (J X ), β, & γ are the values of σv , (J X ), β, & γ that jointly maximize the likelihood at a given m χ .
(J X ), β, and γ are the value of (J X ), β, and γ that jointly maximize the likelihood for a given m χ and σv . We should note that the maximum likelihood values for β and γ are consistent with what we would expect from galactic foreground model of the γ-rays from cosmic-rays [92,93]. The 1D confidence intervals in σv at the nσ confidence level are determined by identifying the range of σv such that ∆ ln L(m χ , σv ) ≤ n 2 /2 .
In the next section we will present the upper limit of the 2σ confidence intervals (95.4% confidence level) for our 3 different selection cuts in the number of clusters. In our analysis, we did not find any significant deviations from background. We found a few regions with 2 σ fluctuations (i.e. consistent with purely background). Figure 10 illustrates out main results. We plot 95% confidence upper limits on dark matter cross section as a function of particle mass for χχ → γγ. The lowest three curves illustrate the results of our likelihood analysis for our 3 cases: on all the clusters, the 2500 brightest clusters, and the 600 brightest clusters. For comparison, we also show the Anderson et al. upper limits [68] (prior to our work, the strongest bounds from clusters). Finally, the figure also includes upper confidence limits on the same clusters as Anderson et al. using our stacking method. Our results are stronger than previous cluster based results, using [68] as a recent comparison point, but not as strong as previous galactic center based dark matter constraint results [94] [95]. As expected from Eq.(2), the bounds on the cross section get weaker for higher WIMP mass.

VI. RESULTS
We note the two most significant upward fluctuations are at around 90 GeV and a bump at 350 GeV, both with a significance of roughly 2 σ. The fluctuation is slightly more pronounced in the case where all the clusters are included, compared to the sub-sample constraints from the 2500 and 600 cluster sets. Just as a matter of total speculation, we remind the reader that 350 GeV is around half the energy of the anomalous bump found at LHC.
A few notes on the relative merits of the different cuts on the clusters. For most values of DM mass, the constraints are stronger when including more galaxy clusters. However, at some masses, the constraints are better for smaller subsets of galaxy clusters. We attribute this to statistical fluctuation. For instance at m χ ∼ 200 GeV, the 600 halo case has stronger constraints than the all halo case. For the all halo case at this energy, there are ∼70 photons (which we assume are background). Given this number, for the 600 halo case, we would expect around 12 photons, but there are only 4 photons. For this case we find that the signal-to-noise ratio is 1.4 times larger for the 600 halo case and thus the bound is stronger.
The top two curves of Figure 10 contrast our stack method with the prior analysis of Anderson etal [68], which studied only 16 clusters. The top two curves in the figure illustrate bounds on the same 16 clusters using our technique (green curve) and the Anderson etal technique (blue curve). Their bounds are expected to be stronger because they obtain a separate likelihood for each of the clusters individually whereas we stack the clusters. In particular, an individual likelihood for each cluster can incorporate spatial information to reduce the background from point sources and other backgrounds, whereas in our stacking method we must instead throw out clusters coincident with the 150 brightest Fermi-LAT γ-ray point sources. Most of the point source photon energies are below 100 GeV, and one can see that the biggest difference between the two approaches (green and blue curves) is at these lower energies. Indeed our results are weaker by up to an order of magnitude, yet do roughly follow the same shape as a function of WIMP mass. In principle the method of [68] does provide stronger bounds, and in the future we plan to use that approach in studying a few hundred clusters from the Tully catalog.
Our final results (the bottom three curves in the figure) provide far stronger bounds than the previous work of [68] due to the simple fact that we are studying 26,000 clusters (rather than the 16 of [68]).

VII. SUMMARY
In this paper we used γ-ray observations of some 26,000 galaxy clusters in Pass-7 Fermi Large Area Telescope (LAT) data to place some of the strongest bounds to date on the cross section of dark matter particles for annihilation to a gamma-ray line, χχ → γγ. The clusters were selected from the Tully 2MASS Groups catalog. For each cluster, we defined a region of interest containing 95% of the expected DM annihilation luminosity. We slid a bin of energy range ∼ twice the energy resolution of Fermi-LAT across our full spectrum of interest. We then added up all the observed photons for all the clusters for each energy window, and compared our observed photon count to the expected count. We searched for a bump in the observed photon count above an expected power law background, i.e. a line (or internal bremsstrahlung [70,76]) signal at an energy equal to the WIMP mass due to DM annihilation. Since no excess above background was found, we used the null signal to place bounds on the DM annihilation cross section as a function of WIMP mass for m χ between 20 and 500 GeV. We found 2σ upward fluctuations most prominently at ∼350 GeV, but these are clearly not significant enough to claim detection.
We have improved on the previous best limits provided from galaxy clusters by a factor of 5 to 10 (depending on the DM mass). Our cluster based constraints are not yet as strong as bounds placed using the Galactic Center, although a less conservative "boost factor" from cluster substructure than the one we have chosen could strengthen our bounds considerably. Our analysis, given this choice of possible boost, is not yet sensitive enough to fully rule out typical realistic DM candidates, especially if the gamma-ray line is not a dominant annihilation mode.
In the analysis in this paper we stacked the clusters. In the future, a stronger bound may be obtained by performing a different analysis, namely determining an individual likelihood function for each of several hundred clusters (this approach was previously used by [68] for 16 clusters). The latter technique can incorporate spatial information to reduce the background.
Additional work for the future would be to repeat our analysis with Fermi-LAT Pass 8 data instead of Pass 7 data, which may yield stronger constraints (the Pass 8 data were not available when we started this work). One could also perform a joint likelihood analysis of both the nearby clusters in conjunction with the Galactic Center.
Our improvement on the previous limits provided from galaxy clusters by a factor of 5 to 10 shows promising future use of galaxy clusters as a viable source for placing further bounds on the particle physics of dark matter. On general merits, by looking at clusters of galaxies which have been characterized comprehensively by Tully, we now have the strongest constraints on dark matter annihilation to a line in clusters of galaxies. In future work, we also plan a similar study on other DM annihilation channels, which will produce a broader range of photon energies.