What the Milky Way's Dwarfs tell us about the Galactic Center extended excess
Ryan E. Keeley, Kevork N. Abazajian, Anna Kwa, Nicholas L. Rodd, Benjamin R. Safdi
UUCI-TR 2017-11MCTP 17-19, MIT-CTP 4943
What the Milky Way’s Dwarfs tell us about the Galactic Center extended excess
Ryan E. Keeley, ∗ Kevork N. Abazajian, † Anna Kwa, ‡ Nicholas L. Rodd, § and Benjamin R. Safdi ¶ Department of Physics and Astronomy, University of California, Irvine, Irvine, CA 92697 Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139 Michigan Center for Theoretical Physics, Department of Physics, University of Michigan, Ann Arbor, MI 48109
The Milky Way’s Galactic Center harbors a gamma-ray excess that is a candidate signal ofannihilating dark matter. Dwarf galaxies remain predominantly dark in their expected commensu-rate emission. In this work we quantify the degree of consistency between these two observationsthrough a joint likelihood analysis. In doing so we incorporate Milky Way dark matter halo profileuncertainties, as well as an accounting of diffuse gamma-ray emission uncertainties in dark matterannihilation models for the Galactic Center Extended gamma-ray excess (GCE) detected by the
Fermi Gamma-Ray Space Telescope . The preferred range of annihilation rates and masses expandswhen including these unknowns. Even so, using two recent determinations of the Milky Way halo’slocal density leave the GCE preferred region of single-channel dark matter annihilation models tobe in strong tension with annihilation searches in combined dwarf galaxy analyses. A third, higherMilky Way density determination, alleviates this tension. Our joint likelihood analysis allows us toquantify this inconsistency. We provide a set of tools for testing dark matter annihilation models’consistency within this combined dataset. As an example, we test a representative inverse Comptonsourced self-interacting dark matter model, which is consistent with both the GCE and dwarfs.
PACS numbers: 95.35.+d,95.55.Ka,95.85.Pw,97.60.Gb
I. INTRODUCTION
The Large Area Telescope aboard the
Fermi Gamma-Ray Space Telescope , Fermi -LAT, has observed a brightexcess of gamma rays towards the Galactic Center whosepresence is robust to systematic uncertainties in the stan-dard background templates [1–14]. This excess has gen-erated a great deal of interest since dark matter (DM)annihilation models can explain three compelling coinci-dences in the signal. First, the excess’ spatial morphol-ogy matches the predictions of a generalized Navarro-Frenk-White (NFW) profile, which is a generic predic-tion of cold DM models [15, 16]. Second, the total fluxof the signal is well fit by the annihilation cross-sectionrequired by a thermal production scenario to generatethe observed cosmological relic abundance. Third, thespectrum roughly matches the expectations of a tens ofGeV weakly interacting massive particle (WIMP) annihi-lating to standard model particles. Should the GCE turnout to be explained by such an annihilating WIMP DMparticle, it could be the first non-gravitational evidenceof DM and the first strong clue of the particle nature ofDM.The prompt annihilation of WIMPs is not the onlyclass of DM models that can explain the GCE, however.For example, a class of self-interacting DM (SIDM) mod-els can explain the GCE via up-scattering of starlightthat would not be seen in dwarf galaxies [17–21]. Specif- ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] ically, this class of SIDM particles could annihilate intoelectrons (as well as the other standard model leptons)and these electrons could up-scatter the Galactic Cen-ter’s interstellar radiation field (ISRF) via the inverseCompton (IC) process.There are also reasonable astrophysical interpretationsof the GCE. Most notably is that the GCE can arise froma population of unresolved millisecond pulsars (MSP)[5, 8, 22–31]. Specifically, observations of MSPs in glob-ular clusters show they have a spectrum consistent withthe spectrum of the GCE. Further, low mass X-ray bina-ries (likely progenitors of MSPs) in M31 have been ob-served to follow a power law radial spatial distribution,similar to the expectations of an NFW halo [5, 27]. Otherastrophysical explanations might include more dynamicevents such as cosmic-ray injection into the Galactic Cen-ter (GC) [32–35]. Furthermore, the presence of the Fermi
Bubbles tell us that such dynamic events have occurredin the past, so whatever mechanism produced the
Fermi
Bubbles, could also have produced the GCE [36–38].There have arisen a number of independent avenuesthat each has the potential to challenge a DM inter-pretation of the GCE. One such avenue is to look fora gamma-ray excess from other DM halos. Such halosinclude those of galaxy clusters, the limits from whichhave recently been extended to be in slight tension withthe GCE [39, 40], and the Milky Way’s satellite dwarfgalaxies, the most recent
Fermi limits appearing in Ref.[41]. Unfortunately,
Fermi -LAT observations of both ofthese sources, and particularly the dwarfs, have not seena significant complementary gamma-ray excess. In par-ticular, this difference between the GCE and the dwarfs Note, however, there has been a low-significance detection of a a r X i v : . [ a s t r o - ph . H E ] O c t has the potential to rule out certain classes of DM inter-pretations of the GCE. Specifically, any minimal modelbased around a two-body annihilation process (any pro-cess where the flux is proportional to the square of theDM density) would exhibit this same tension.Other avenues to test whether the GCE is better ex-plained by annihilating DM or astrophysics is to moreprecisely check whether the morphology of the excesstruly follows a smooth NFW profile. Tension in themorphology of the GCE has arisen from the detection ofan ‘X-shape’ residual in the Fermi data which correlatesinfrared emission as seen by the WISE telescope [45].Should this ‘X-shape’ template account for the entirety ofthe GCE, it would challenge any DM interpretation sinceDM annihilation would not produce such a shape. Mea-surements of the GCE being consistent with wavelets [46]or non-Poissonian fluctuations have also been reported[47–50], which would indicate an MSP rather than DMorigin, though systematic uncertainties in such analysesremain [51]. Specifically, small scale gas clouds are leftout of the model used by
GALPROP , software used to gen-erate the gamma-ray templates associated with cosmicrays propagating through the Milky Way, which couldconfuse any detection of point sources near the GC [51].Though each of these lines of evidence against a DM in-terpretation of the GCE have their own systematic un-certainties, many of these systematics are independentof each other. Arguably, these different lines of evidenceadd up to strongly indicate that the GCE is astrophysicalin origin.Our focus in the present paper is to consider one as-pect of this general line of reasoning: the consistency be-tween the GCE and the dwarfs, and to do so with a moredetailed treatment of the systematics coming from bothsides. The discussion is structured as follows. In sectionII we discuss the background models we investigate tounderstand some of the dominant sources of systematicuncertainties in the problem. In section III, we discussDM annihilation models of the GCE. In section IV, wediscuss alternative models to promptly annihilating DM,including astrophysical interpretations and SIDM mod-els. We conclude in section V.
II. BACKGROUND MODELS AND DATA
There remains significant uncertainty regarding thevarious processes that contribute to the gamma-ray sig-nal coming from the GC. Any interpretation of the GCEwill necessarily be affected by these uncertainties. Tocapture these effects, we investigate four different cases ofthe astrophysical background contributions to the GCE.For all of our cases (denoted cases A, B, C, and D),we use data collected by
Fermi -LAT. For cases A, B, gamma-ray excess from Reticulum II and Triangulum II, see e.g.[42, 43], although see also [44]. and C, that data corresponds to observations over a 103month period from August 2008 to March 2017. We useall
SOURCE -class photons from the Pass 8 instrument re-sponse functions. We apply a maximum zenith angle cutof 90 ◦ to avoid contamination. In cases A, B, and C,we focus our analysis on the innermost 7 ◦ × ◦ region ofinterest (ROI) about the Galactic Center. We then binthese photons into spatial bins of size 0 . ◦ × . ◦ foreach energy bin. The photon events range from 200 MeVto 50 GeV and are binned in 16 logarithmically spacedenergy bins.For Case D we instead choose a dataset similar to thatconsidered in the Inner Galaxy analyses of [10, 49]. Herewe use the best quartile, as graded by the Fermi pointspread function, of
ULTRACLEANVETO -class Pass 8 pho-tons, gathered between August 4, 2008 and June 3, 2016with recommended quality cuts. This case can be con-trasted with the above in that it generally correspondsto less data, but with less cosmic-ray contamination andimproved angular reconstruction per event. To mimican earlier Inner Galaxy analyses, we use a larger ROIof 30 ◦ × ◦ , masking latitudes less than 1 ◦ . No maskswere applied to the data in cases A-C. We also maskthe top 300 brightest and most variables sources in the3FGL catalog [52] at 95% containment. The photons arebinned into 40 equally spaced logarithmic bins between200 MeV and 2 TeV, and spatially using an nside =128 HEALPix grid [53].With this processed data, we perform a maximumlikelihood analysis to determine the best fit backgroundmodel and GCE model. For each component of ourmodel, we generate a template which encodes the spa-tial distribution of the photons for that component. Thequantity that we are trying to determine is then thelinear combination of these spatial templates that bestfit the observed number of counts. The templates fallinto three groups: point sources, extended emission, anddiffuse emission. The point sources we use are takenfrom
Fermi ’s 3FGL point source catalog [52] and theyare typically well characterized or independent of theGCE result. The extended emission components includea GCE template, as well as background components com-ing from cosmic rays interacting with gas in the interstel-lar medium (ISM) or photons in the ISRF. Specifically,these would include any IC radiation from high energyelectron cosmic rays up-scattering the ISRF, neutral pion( π ) decay generated from hadronic cosmic rays interact-ing with the ISM, and bremsstrahlung radiation arisingfrom high energy electrons interacting with the ISM. Thespatial distribution of these components are more a priori uncertain than point sources and are partially degeneratewith the GCE, especially in the lowest energy bins, wherethe point-spread function is the largest. Therefore, it isthese uncertainties and degeneracies that make a carefuland broad investigation of the diffuse backgrounds cru-cial to analyzing the GCE and are the main differencebetween our different cases.Since the uncertainty in the spectral shape of the GCE Energy [GeV] -8 -7 -6 -5 E d N / d E [ G e V c m − s − ] Case A
GCE π + bremsICTotal ModelCounts Energy [GeV] -8 -7 -6 -5 E d N / d E [ G e V c m − s − ] Case B
GCEp8v620 cm brems3.4 µ m ICTotal ModelCounts Energy [GeV] -8 -7 -6 -5 E d N / d E [ G e V c m − s − ] Case C
GCEp7v620 cm brems3.4 µ m ICTotal ModelCounts Energy [GeV] -8 -7 -6 -5 E d N / d E [ G e V c m − s − ] Case D
GCEp6v11IsotropicBubblesTotal ModelCounts
FIG. 1. Here we plot the energy flux spectrum (intensity) E dN/dE for the various templates included in the likelihood fitsfor our A, B, C, and D background cases. These show the total emission from the ROI, 7 ◦ × ◦ for cases A-C and 30 ◦ × ◦ for case D. The error bars on the counts is the Poisson error. The various 3FGL sources were also varied in the fits but are notincluded for the sake of simplicity. signal is dominated by systematic uncertainties in thebackground templates, rather than Poisson fluctuationsof the total counts, it is necessary to explore multiplepossible background models. To this end, we use fourdifferent sets of templates for these extended backgroundmodels: • Case A:
We use the templates for the π ,bremsstrahlung, and IC emission for ‘model F’from Horiuchi et al. (2016) [51], which in turn used diffusion model parameters from Calore et al.(2014) [11] to generate their background models.Their ‘model F’ corresponds to the diffuse back-ground model that was found to best fit the datain their ROI. Unlike the Fermi collaboration Pass8 and Pass 7 diffuse backgrounds, the IC compo-nent of the diffuse background is fit independentlyof the π +bremsstrahlung components. We usedthe templates for ‘model F’ from these papers. Inthis work, we denote this ‘case A.’ • Case B : For this case, we use the Pass 8 Galacticinterstellar emission model from the
Fermi tools,which models the distribution of gamma rays com-ing from π decay, bremsstrahlung, and IC. Allthree components are combined in a single diffusetemplate with fixed relative normalizations in eachenergy bin. Furthermore, we used a template whichtraces the 20 cm radio emission first discovered byYusef-Zadeh et al. (2013) [54]. We also include atemplate for an additional IC component that wasderived from 3.4 µ m maps from the WISE telescope,discovered by Abazajian et al. (2014) [12]. • Case C:
This case uses the same templates for thebremsstrahlung and IC components but uses p7v6 model for Pass 7. • Case D:
This case uses the p6v11 template andfloats an isotropic template as well as a templatefor the
Fermi bubbles.For all these cases, we allow the flux associated witheach template in a given bin to be independent of the fluxin other bins, rather than assume a specific componenthas a specific spectral shape. This allows us to be ag-nostic about the shape of the spectrum for each of thesesources, but potentially comes at the cost of over-fittingthe data. The results of these maximum likelihood fitsfor cases A-D are shown in Fig. 1.To calculate posteriors for the dwarfs, we use the fluxlikelihood limits from Albert et al. (2016) [41]. Specifi-cally, we use the flux likelihood manifolds for the nineteenkinetically confirmed dwarf galaxies that have measuredJ-factors. The J-factor for Reticulum II is calculated inSimon et al. (2015) [55] and the rest are calculated inGeringer-Sameth e t al. (2014) [56]. To calculate theseflux likelihood limits, Albert et al. use six years of LATdata with 24 equally-spaced logarithmic bins between500 MeV and 500 GeV. They binned the photons in a10 ◦ × ◦ region about the target dwarf galaxies with apixel size of 0 . ◦ in order to model any overlap from thepoints spread function of the point sources in the 3FGLcatalog, from the Galactic diffuse emission, and from theisotropic model. Each target dwarf galaxy was modeledas a point-like source and used a maximum likelihoodanalysis with these templates to generate the flux likeli-hood limits. III. ANNIHILATING DARK MATTER MODELSA. Flux Spectra
The differential flux in some ROI for the class of two-body DM annihilation is given by the following: d Φ dE = 14 π Jm χ (cid:104) σv (cid:105) dNdE . (1) Here, J is the J-factor, the integral of the density-squared over the ROI and through the line of sight. m χ is the mass of the DM particle and (cid:104) σv (cid:105) is the ther-mally averaged cross section of the annihilation. dNdE isthe per-annihilation spectrum, which we calculated using PPPC4DMID [57]. For our dark matter models, we use flatpriors on the DM mass and scale invariant priors on thecross section. The prior on the J-factor is discussed inthe next section.
B. J-factors
The J-factor is the square of the DM density integratedthrough the line of sight and integrated over the ROI. J = (cid:90) ROI d Ω (cid:90) dz ρ ( r ( z, Ω)) . (2)As in Abazajian & Keeley 2015 [58], we determine theprior on the J-factor for the GC by parameterizing theMilky Way’s DM halo as a generalized NFW profile witha local DM density ( ρ (cid:12) ), a scale radius ( R s ), and an innerprofile slope ( γ ) ρ ( r ) = ρ (cid:12) (cid:16) rR (cid:12) (cid:17) γ (cid:16) r/R s R (cid:12) /R s (cid:17) − γ . (3)Each of these parameters has a probability distribution,so in principle, we could say the prior on the J-factoris the product of the probability distributions of eachof these parameters and then perform the change of vari-ables to write this probability distribution as a function ofthe J-factor. This is analytically cumbersome, so we usenumerical Monte Carlo methods to calculate this distri-bution. Specifically, we draw values for the local density,scale radius, and inner slope to compute a set of J-factorsand then use kernel density estimation to define the priorfor the GCE J-factor.For the local density, we use the value determined byZhang et al. (2012) [59]: ρ (cid:12) = 0 . ± .
08 GeV cm − .This robust determination of the local DM density is de-rived from modeling the spatial and velocity distributionsfor a sample of 9000 K-Dwarf stars from the Sloan Digi-tal Sky Survey (SDSS). The velocity distribution of thesestars directly measures the local gravitational potentialand, when combined with stellar density constraints, pro-vides a measure of the local DM density.The prior on the scale radius is calculated from the con-centration, which is the ratio of the virial radius to thescale radius. The uncertainty in the concentration is cal-culated from simulations of galaxy formation. Sanchez-Conde and Prada (2013) [60] parameterized the uncer-tainty in the concentration of a DM halo as a functionof that halo’s mass. Thus we can write the prior on thescale radius as:log L = − (log ( R vir /R s ) − log c ( M vir )) × . . (4) J-factor [GeV cm − sr]0.00.20.40.60.81.0 S ca l e d P r ob a b ilit y Case ACase BCase CCase D
FIG. 2. The prior on the J-factor integrated over the ROIderived through a Monte Carlo convolution of the priors onthe local density, scale radius, and inner slope. Since each ofthe different background cases have different best fit values for γ , and since case D corresponds to a larger ROI, the deriveduncertainties on the J-factors are different. The prior on the inner slope we use for the MonteCarlo calculation of the J-factor is taken to be the pos-terior determined by the spatial information containedin the GCE data. We constrain the inner slope by run-ning the likelihood analysis with the same backgroundmodels but with different NFW spatial templates thathave different values for the inner slope. The likelihoodanalysis calculates the ∆Log-likelihood value for each ofthese different cases, which allows us to fit a χ distribu-tion to these ∆Log-likelihood values. This determines thebest fit value of γ and its error. Unsurprisingly, this de-rived prior on the inner slope depends on the backgroundmodel used. For case A we calculate, γ = 1 . ± .
04; forcase B, γ = 1 . ± .
04; for case C, γ = 1 . ± .
04; andfor case D, we calculate γ = 1 . ± .
06. The results ofthis Monte Carlo calculation of the priors on the J-factorfor the different background cases is shown in Fig. 2.We employ the priors on the J-factors for the dwarfgalaxies from Albert et al. (2016) [41]. These are allreported as log-normal distributions. These J-factorscome with some caveats, however. Specifically, assump-tions about how spherically symmetric the dwarf galaxyis, which in turn can influence the inferred cuspiness ofthe density profile, can lead to systematic uncertaintiesgreater than the statistical uncertainties [61–63].
C. Evidence Ratios
To quantify the tension between the GCE and thedwarfs, we calculate a Bayesian evidence ratio. This ev-idence ratio can be used in answering the question: bywhat factor do the odds of some model being true changewith the inclusion of a new data set. It is the product of the Bayesian evidences of two data sets, D and D ,when considered separately divided by the evidence ofthe two data sets when considered jointly [64]:ER = p ( D ) p ( D ) p ( D , D )= (cid:82) dθ p ( D | θ ) p ( θ ) (cid:82) dθ p ( D | θ ) p ( θ ) (cid:82) p ( D , D | θ ) p ( θ ) . (5)This can be interpreted as a Bayes factor where the twomodels being compared are the same except for the factthat the model corresponding to the numerator has anadditional, independent copy of the parameter space andthe two parameter spaces describe the two data sets sep-arately.This statistic can indicate three different outcomes forthe model. First, if the data set D contains no informa-tion, then this evidence ratio is unity. If D is entirelyconsistent with D then the evidence ratio should be lessthan unity. This is expected since increasing the com-plexity of a model should come at a cost of subjectivebelief. If D is in tension with D then this evidence ra-tio will be greater than unity. How strongly this evidenceratio prefers consistency or tension can be interpreted byany standard Bayes factor scale. In this work, we chooseto interpret our evidence ratios by the Jeffreys’ scale.Using this setup, we then calculate evidence ratios be-tween the combined dwarf galaxies and the GC. The re-sults are stated in Table I.One particularly useful feature of evidence ratios inthis context is that, compared to Bayes factors, they arerelatively insensitive to systematic uncertainties in thebackground models. These systematic uncertainties canalter the total flux of the signal, but they more drasti-cally change in which energy bin this flux is distributed.This is seen most clearly in the lowest energy bins, wherethe inclusion of diffuse templates from 20 cm maps ofbremsstrahlung emission and 3.4 µ m maps of IC emis-sion, for cases B and C, remove all the photons from theNFW template for these bins. Such changes to the low-est energy bins changes the overall curvature of the GCEspectrum, which, in turn, significantly changes the bestfit mass but not the best fit cross section [58]. Whenthe best fit mass of the GCE changes, the amount ofoverlap between the GCE posterior with the combineddwarfs posterior (and hence the evidence ratio) changesrelatively little. This lack of change in overlap comesfrom the fact that the contours of the dwarf posteriorare almost parallel to contours of constant cross section,since the lack of a dwarf signal contains no significantamount of information about the spectrum. It is becausethe evidence ratio is most sensitive to the cross sectionand not particularly sensitive to the dark matter massthat the evidence ratios are more robust to systematicuncertainties in the background templates. This is bornout in Table I where the DM evidence ratios for the dif-ferent cases vary by only two orders of magnitude. Onthe other hand, because the Bayes factors are sensitive TABLE I. Evidence ratios for our five models using the diffusetemplates for our various background cases.Model Case A Case B Case C Case DDM: b ¯ b τ + τ − to both the normalization and the shape of the spectrumit can vary by 30 orders of magnitude, as seen in TableII.Beyond systematic uncertainties due to the inclusionof additional templates for bremsstrahlung and inverseCompton processes, uncertainties in the diffuse model forthe GALPROP generated π , IC, and bremsstrahlung tem-plates, can alter the total flux of the GCE signal and willaffect the best fit cross section for the GCE and henceaffect the tension with the dwarfs. This is seen by thefact the evidence ratio for our different cases significantlychange. This change is caused by the fact that the differ-ences in these diffuse emission templates, for cases B andD, shift the overall flux of the GCE signal to smaller val-ues, relative to case A. In these background model cases,the presence of the GCE is less significant and reducesthe significance of the tension with the dwarfs.The information encoded by the evidence ratio can bequalitatively seen in Fig. 3, which plots the posterior ofour b ¯ b and τ + τ − DM annihilation models for each of ourdifferent GCE background cases and for the dwarf data.The amount of overlap in the GCE posteriors and thedwarf posterior indicates the amount of tension betweenthe data sets.For our DM annihilation models, we calculate evidenceratios between 15 and 2200 for the b ¯ b channel and be-tween 27 and 4300 for the τ + τ − channel. Using theJeffreys Scale, this indicates a strong ( >
10) to decisive( > b ¯ b or τ + τ − . Any model of prompt two-body decay,described by a J-factor, would exhibit this same tension.Hence, models that contain only novel versions of spec-trum dN/dE , or branching ratios, will not alleviate thisstrong tension. D. Caveats
The GCE-dwarf tension we quantified in the previ-ous section certainly depends on the prior informationadopted for the J-factors of the GC region and the dwarfgalaxies. Naturally, if there was a significant change in the inferred DM content of either the GC region or dwarfgalaxies, then the nature of tension would correspond-ingly change. However, our choices for the J-factor ofthe GC region and dwarf galaxies are those determinedby the most robust analyses available.The parameter that the J-factor is most sensitive tois the local density of DM. As stated in a previous sec-tion, we use a value of 0.28 ± taken fromZhang et al. (2012) [59]. Other groups including Patoet al. (2015) [65] and McKee et al. (2015) [66] tend tofind higher values for the local density. To fully resolvethe tension between the GCE and the dwarfs, the GCEJ-factor needs to increase between 1 and 1.5 orders ofmagnitude, which translates into a local density of 3 to 6times greater. As we show, none of these determinationsof the local density relieve the GCE-dwarf evidence ratioto be unity.Another parameter with a systematic uncertainty isthe scale radius of the DM profile. Small deviationsaround our fiducial value would not change the J-factorby a great deal since the inner profile is unchanged dueto the scale radius being beyond the local radius. How-ever, should the scale radius become smaller than thelocal radius, the inner density profile would increase as r − between the local radius and the scale radius, result-ing in a larger J-factor. A profile with such a small scaleradius could only occur in halos with a concentration pa-rameter far outside of what CDM simulations predict forhalos with the mass of the Milky Way.The inner slope γ is more robust to systematic uncer-tainties, in that it is determined directly from the spatialinformation of the gamma-ray data. In particular, tofully resolve the GCE-dwarf tension, a value of around γ = 1 . γ = 1 . .
3. Despite systematic uncertainties in theparameterization of the Milky Way’s DM profile, no sin-gle alteration can fully relieve the tension between theGCE and dwarf data.
IV. MODELS
We have shown that there is tension with the standardWIMP scenario between the derived cross sections fromthe GC and the dwarfs, with some important caveats.This tension can potentially point to alternate modelsbeing better explanations for the GCE, including astro-physical interpretations to more complicated DM mod-els. To quantitatively answer this question, we calculatea Bayes factor: K = p ( M | D ) p ( M ) p ( M | D ) P ( M ) = p ( D | M ) p ( D | M ) . (6)We consider the following models: two astrophysical in-terpretations, one with a log-parabola spectrum and an-other with an exponential cutoff spectrum, and a SIDM
30 40 50 60 70 80 90 100Mass [GeV]10 -27 -26 -25 -24 -23 C r o ss S ec ti on [ c m s ec − ] b ¯ b DwarfsABCD -27 -26 -25 -24 -23 C r o ss S ec ti on [ c m s ec − ] τ + τ − DwarfsABCD
FIG. 3. Here we show the 1, 2, and 3 σ contours of the posteriors for the annihilation cross section and DM mass. Our calculatedlimits on the dwarf signal is in green, case A is in orange, case B is in blue, and case C is in pink. The results for b ¯ b on the leftand τ + τ − on the right. The amount of overlap qualitatively demonstrates the information contained in the evidence ratio andshows how consistent two-body DM annihilation models are at explaining both the GCE and the lack of a dwarf signal. model where the GCE gamma rays are generated by DMdecaying to high-energy electrons up-scattering starlight.The Bayes factors for our models are given in Table II. A. Astrophysical Interpretations
Should the GCE have an astrophysical interpretation,the gamma-ray spectrum can be parameterized as a log-parabola or a power law with an exponential cutoff. Weinvestigate both parameterization as explanations of theGCE.The spectrum for our log-parabola model is given by: dNdE = N (cid:18) EE s (cid:19) − α − β log( E/E s ) , (7)where N is an arbitrary normalization, E s is a scale en-ergy, α is the slope of the power-law part of the spectrum,and β parameterizes the turnover of the spectrum.The spectrum for our power law with an exponentialcutoff model is given by: dNdE = N (cid:18) EE s (cid:19) γ e − E/E c , (8)where N is the normalization of the spectrum, E s is ascale energy, γ is the slope of the power-law part of thespectrum, and E c parameterizes how fast the spectrumcuts off.Our astrophysical models do not have a specific phys-ical interpretation so it is not straightforward to inves-tigate to what extent the GCE and the lack of signalfrom the dwarf galaxies are compatible given these mod-els. Presumably, if the GCE and any potential dwarf signal were to be explained by the same category of astro-physical object, then they should have the same spectralparameters. Therefore, it makes sense for our model tohave only one set of spectral parameters that describesboth the GCE and the dwarfs. The normalizations of thespectrum, however, would not necessarily be the same.One option is to allow the normalization of the spectrumof the GCE and the spectrum of each of the dwarfs to beindependent. Following this parameterization, we calcu-late an evidence ratio between the GCE and the dwarfsof about 1, which would indicate the two data sets con-tain no new information relative to each other. This isexpected, since if we put in the fact that the signals areindependent, we should get out that they have no mutualinformation. Instead of saying these normalizations areentirely independent of each other, we use a zeroth orderansatz to parameterize the normalization as the productof the stellar mass of the system and the gamma-ray rateper stellar mass. The stellar mass would, of course, beindependent between regions, but the gamma-ray rateper stellar mass should be the same between regions. Tothis end, we find N in the above equations such thatthe integral of dN/dE over our energy range (200 MeVto 50 GeV) is one. This allows us to attach physicalinterpretations to our normalization for d Φ /dE .Specifically, it makes sense, should the initial massfunction of some galaxy be independent of the stellarmass of that galaxy, that the gamma-ray production ratescales linearly with the stellar mass of the galaxy. Hence,the gamma-ray rate per stellar mass should be consistentacross all regions.Ultimately, this leads to the following parameterizationof the differential number flux: d Φ dE = ˙ N πR M ∗ M dNdE , (9)where M ∗ is the stellar mass of the object, R is the dis-tance to the object, and ˙ N /M is the gamma-ray rate perstellar mass, which should be the same between differentobjects.For both spectra of astrophysical models, we marginal-ized over the spectral parameters with flat priors, andmarginalized over the over the gamma-ray rate per stel-lar mass with a scale invariant prior. We use the stel-lar mass of the dwarfs, the distance to them, as well asthe uncertainties in those parameters from McConnachie(2012) [67]. Interestingly, both of our astrophysical mod-els pick out values for the gamma-ray rate per stellarmass around 10 ± s − M − (cid:12) , which is consistent withknown millisecond pulsars. In the end, the evidence ra-tios for each of our two spectral choices for astrophysicalmodels, for all of our background cases, are less thanunity. Importantly, this less than unity evidence ratioindicates that the combined dwarf and GCE data havea weak indication of a mutual astrophysical excess de-scribed by a single set of parameters.The Bayes factors we compute also point towards apreference for these astrophysical models. As seen in Ta-ble II, the log-parabola spectrum is preferred over anyDM model in each of the cases, and the exponential cut-off spectrum is preferred in three out of four of the cases.The preference in the Bayes factor can be thought of ascoming from two distinct sources. One is the GCE dataon their own prefer that model and the other is thatthe model can better explain the differences in the fluxfrom the GC and the dwarfs. Astrophysical interpreta-tions, with evidence ratios less than unity, can do betteron the latter count, but interestingly, depending on thedata case, can also do better on the former count. In allcases, the log-parabola spectrum can explain the GCEdata better than dark matter models, but in cases B andC, the exponential cutoff can do so also. This preferencein some cases for the log-parabola spectrum is predom-inantly coming from the lowest energy bins. The maxi-mum likelihood fit prefers giving no appreciable amountof photons to the lowest energy bins, a fact that is diffi-cult for DM models to explain, but is more easily accom-modated by the log-parabola spectra. This preference ofthe lowest energy bins for the log-parabola spectra canbe seen in Figure 4, where we plot the best fit models,along with the data. It is worth noting that these lowestenergy bins have the largest systematics associated withthem due to the large size of the size of the point spreadfunction at those energies [11]. Unlike the evidence ra-tios, the Bayes factors are particularly sensitive to thesesystematics, particularly because no model is a strikinglygood fit, just less bad than the others. Indeed, when ig-noring the first few data points for each data case, theBayes factors tend to show less extreme results, givingmore consistent fits to the GCE. With these truncateddata sets, the values of the Bayes factors come from themodels’ abilities to explain the difference in flux comingfrom the GC and the dwarfs.On an additional note, the preference for b ¯ b can also be seen in Fig. 4. Since the τ + τ − model requires a light(compared to the b ¯ b model) dark matter mass to explainthe peak of the GCE spectra at around 1-2 GeV, andsince these annihilating dark matter models cutoff in en-ergy at around their mass, the τ + τ − model fail to ac-count for the GCE spectra that gradually fall off withlarge energies, such as in cases A, B, and D. B. A Representative SIDM Model
In certain classes of SIDM models for the GCE, thegamma-ray excess is generated by the DM particles an-nihilating to electron-positron pairs through a light me-diator [17]. The electrons then up scatter the starlight inthe Galactic Center via an IC process. This would nat-urally explain the difference in the observed gamma-rayflux between the GC and the dwarfs since the stellar den-sity of the dwarfs, and therefore the interstellar light, ismany orders of magnitude smaller than the stellar massof the GC.Should the GCE be explained DM annihilating to elec-trons that interact with the ambient starlight, the processshould be governed by the following IC equation: dn γ dEdt = σ T cn e n ISRF dN γ dE , (10)where n γ is the number density of gamma rays, σ T is theThomson cross section, n e is the number density of elec-trons produced by annihilating DM, n ISRF is the numberdensity of low energy photons in the interstellar radiationfield, and dN γ dE is the probability distribution function ofproducing a gamma ray of energy E via this IC process.Naturally, this probability distribution function dependson the probability distribution functions of the energiesof the electrons produced via DM annihilation and theenergy distribution of the starlight: dN γ dE = (cid:90) dE e dE ISRF p ( E γ | E e , E ISRF ) p ( E e ) p ( E ISRF ) . (11)In principle, other energy-loss mechanisms, such as syn-chrotron emission, can alter the energy distribution ofelectrons in this IC process. We checked this modelagainst the spectrum PPPC4DMID calculates and found theshape of the spectra were largely consistent.Since the electrons are produced via a two-body in-teraction, the number density of electrons should scaleas the square of the number density of DM particles: n e ∝ n χ . To convert the time derivative of the differen-tial number density of photons to some number flux, weneed to evaluate the following integral: d Φ γ dE = (cid:90) dV (cid:48) π ( (cid:126)R − (cid:126)R (cid:48) ) dn γ dEdt ( (cid:126)R (cid:48) ) . (12)Choosing the origin of the coordinate system to be at R = 0 leads to the standard expression for the J-factor,should the process be entirely a two-body process andthe time derivative of n γ scale solely as the square of theDM particles. Putting this all together, we get: d Φ γ dE ∝ (cid:90) d Ω dz π n ISRF n χ dNdE . (13)Instead of using this equation as written, we make thefollowing assumptions and simplifications. First, n ISRF is approximately constant where the density of DM islargest, so we can pull the factor of n ISRF outside theintegral. Second, it should be true that the number den-sity of photons from stars scales with the stellar mass ofthose stars we replace n ISRF with the stellar mass of thegamma-ray source, M ∗ : d Φ γ dE ∝ Jm χ M ∗ M ∗ , GC . (14)Taking this spectrum leads to a model that has fargreater consistency between the GCE and dwarfs; theevidence ratios for this model are all around unity foreach of the data cases. This highlights the possibility toalleviate the tension when going beyond simple two-bodyfinal state scenarios.The best fit DM mass for this representative SIDMmodel is 15 ± ± ± V. CONCLUSIONS
We have analyzed the GCE in a wide variety of back-ground models by performing a template based likeli-hood analysis of the GC using four different models forthe diffuse background templates. To answer the ques-tion of whether an annihilating DM interpretation can beconsistent with the lack of dwarf signals, we calculatedevidence ratios for each model of the GCE and for eachcase of diffuse background models. These evidence ratiosare sensitive to the choice of background model but they
TABLE II. Bayes factors for the considered models, relativeto the b ¯ b model, for each of the different background cases.Values larger than one indicate the data prefer that modelover b ¯ b .Model Case A Case B Case C Case DDM: τ + τ − × − × − × × − Log-Parabola 3 × × × × Exponential Cutoff 2 × × × × − × − × − all display strong to decisive tension between the GCEand the dwarfs for annihilating DM models. Specifically,cases A and C show decisive tension, with evidence ra-tios greater than 100 for both annihilation channels, andcases B and D show strong tension with evidence ratiosgreater than 10 for both channels. This difference can,at least in part, be attributed to the fact the likelihoodfit for these cases seem to prefer both giving less flux tothe DM GCE component, and also prefer an NFW tem-plate with a higher value for the inner slope γ . Since thetension is seen to various degrees using a variety of mod-els for the the diffuse emission, it is robust to say thatprompt two-body annihilating DM interpretations of theGCE are in strong doubt.Astrophysical and SIDM interpretations of the GCEfare better with evidence ratios around unity. Ultimately,allowing the gamma ray signal to scale with the stellarmass, as for astrophysical models, or with the productof the J-factor and stellar mass, as with SIDM models,relieves any tension between the GCE signal and lack ofa dwarf signal.We also calculated Bayes factors for our different DMGCE interpretation models. This Bayes factor can bethought of as coming from two different sources: the abil-ity of the model to explain the GCE and the ability of themodel to explain the difference in GCE and dwarf fluxes.These Bayes factors decisively prefer the log-parabolaspectrum model over the DM annihilation models in allof our background cases, and prefer the exponential cut-off model in three of the four background cases. Thispreference for either astrophysical spectrum model pre-dominantly comes from the lowest energy bins where thelikelihood analysis prefers to attribute no amount of fluxto an NFW template. However, these are also the en-ergy bins that have the largest systematic uncertaintiesassociated with them. Standard two-body DM annihila-tion models cannot explain these low energy gamma-raydata, while more general log-parabola and exponentialcutoff models are able to do so. With the long integra-tion time now available from the Fermi -LAT observationsof the GCE, the data allows us to make very precise de-terminations of the GCE’s spectral parameters, given aparticular background model. However, the accuracy of0 Energy [GeV]10 -9 -8 -7 N u m b e r F l ux [ C oun t s c m − s − ] b ¯ bτ + τ − Log-parabolaExponential cutoffSIDMCase A GCE Energy [GeV]10 -9 -8 -7 N u m b e r F l ux [ C oun t s c m − s − ] b ¯ bτ + τ − Log-parabolaExponential cutoffSIDMCase B GCE Energy [GeV]10 -9 -8 -7 N u m b e r F l ux [ C oun t s c m − s − ] b ¯ bτ + τ − Log-parabolaExponential cutoffSIDMCase C GCE Energy [GeV]10 -9 -8 -7 N u m b e r F l ux [ C oun t s c m − s − ] b ¯ bτ + τ − Log-parabolaExponential cutoffSIDMCase D GCE
FIG. 4. Here we plot the number flux for the GCE template along with the best fit spectra for the different models considered.The error bars correspond to the 1- σ region of each bin’s number flux likelihood profiles. these background models are still uncertain. In otherwords, the systematic uncertainties in the backgroundmodel cases dominate over the Poisson statistical uncer-tainties. In fact, there exist two sets of tests. One iswhether DM or astrophysical spectral models’ can ex-plain the joint GCE and dwarf data. The second is theintrinsic ability of the GCE spectral choices to explainthe GC observations. Importantly, the biggest change inthe models’ Bayes factors comes from the spectral mod-els’ different ability to properly fit the GCE. In almostall cases, log-parabola spectra is decisively better in theirevidence ratios at fitting the GCE data (cf. Table II).Therefore, given the GCE data alone, the log-parabolaastrophysical interpretation of the GCE is favored.Furthermore, the combined GCE-dwarf data stronglyto decisively disfavor single channel DM annihilation in-terpretations of the GCE. Secondary-emission from DMmodels like that from SIDM could alleviate the incon-sistent emission between the GCE and dwarf galaxies.Further detailed analysis of the diffuse emission towards the GC will help determine the true nature of the GCEand its relation to any emission from the dwarf galaxies. ACKNOWLEDGMENTS
K.N.A. and R.E.K. are supported in part by NASA
Fermi
GI grant 15-FERMI15-0002. A.K. is supported byNSF GRFP Grant No. DGE-1321846. N.L.R. is sup-ported by the DOE under contracts DESC00012567 andDESC0013999. We thank Dan Hooper, Manoj Kapling-hat and Savvas Koushiappas for useful discussions. Wealso thank Barbara Szczerbinska and the other organizersof
XI International Conference on Interconnections be-tween Particle Physics and Cosmology , in Corpus Christi,Texas, where part of this collaborative work was initi-ated.1 [1] L. Goodenough and D. Hooper, (2009), arXiv:0910.2998[hep-ph].[2] V. Vitale and A. Morselli (Fermi-LAT), in
Fermi gamma-ray space telescope. Proceedings, 2nd Fermi Sympo-sium, Washington, USA, November 2-5, 2009 (2009)arXiv:0912.3828 [astro-ph.HE].[3] D. Hooper and L. Goodenough, Phys.Lett.
B697 , 412(2011), arXiv:1010.2752 [hep-ph].[4] D. Hooper and T. Linden, Phys.Rev.
D84 , 123005(2011), arXiv:1110.0006 [astro-ph.HE].[5] K. N. Abazajian and M. Kaplinghat, Phys.Rev.
D86 ,083511 (2012), arXiv:1207.6047 [astro-ph.HE].[6] D. Hooper and T. R. Slatyer, Physics of the Dark Uni-verse , 118 (2013), arXiv:1302.6589 [astro-ph.HE].[7] C. Gordon and O. Macias, Phys.Rev. D88 , 083521(2013), arXiv:1306.5725 [astro-ph.HE].[8] O. Macias and C. Gordon, Phys. Rev.
D89 , 063515(2014), arXiv:1312.6671 [astro-ph.HE].[9] K. N. Abazajian, N. Canac, S. Horiuchi, and M. Kapling-hat, Phys.Rev.
D90 , 023526 (2014), arXiv:1402.4090[astro-ph.HE].[10] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden,S. K. N. Portillo, N. L. Rodd, and T. R. Slatyer, Phys.Dark Univ. , 1 (2016), arXiv:1402.6703 [astro-ph.HE].[11] F. Calore, I. Cholis, and C. Weniger, JCAP , 038(2015), arXiv:1409.0042 [astro-ph.CO].[12] K. N. Abazajian, N. Canac, S. Horiuchi, M. Kaplinghat,and A. Kwa, JCAP , 013 (2015), arXiv:1410.6168[astro-ph.HE].[13] M. Ajello et al. (Fermi-LAT), Astrophys. J. , 44(2016), arXiv:1511.02938 [astro-ph.HE].[14] C. Karwin, S. Murgia, T. M. P. Tait, T. A. Porter,and P. Tanedo, Phys. Rev. D95 , 103005 (2017),arXiv:1612.05687 [hep-ph].[15] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astro-phys. J. , 563 (1996), arXiv:astro-ph/9508025 [astro-ph].[16] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astro-phys. J. , 493 (1997), arXiv:astro-ph/9611107 [astro-ph].[17] M. Kaplinghat, T. Linden, and H.-B. Yu, Phys. Rev.Lett. , 211303 (2015), arXiv:1501.03507 [hep-ph].[18] M.-Y. Cui, C. Zhang, and H.-S. Zong, Phys. Rev.
D93 ,123516 (2016), arXiv:1605.08138 [hep-ph].[19] T. Lacroix, O. Macias, C. Gordon, P. Panci, C. Bhm, andJ. Silk, Phys. Rev.
D93 , 103004 (2016), arXiv:1512.01846[astro-ph.HE].[20] T. Lacroix, C. Boehm, and J. Silk, Phys. Rev.
D90 ,043508 (2014), arXiv:1403.1987 [astro-ph.HE].[21] I. Galon, A. Kwa, and P. Tanedo, JHEP , 064 (2017),arXiv:1610.08060 [hep-ph].[22] K. N. Abazajian, JCAP , 010 (2011),arXiv:1011.4275 [astro-ph.HE].[23] H. Ploeg, C. Gordon, R. Crocker, and O. Macias, JCAP , 015 (2017), arXiv:1705.00806 [astro-ph.HE].[24] M. Ajello et al. (Fermi-LAT), Submitted to: Astrophys.J. (2017), arXiv:1705.00009 [astro-ph.HE].[25] N. Mirabal, Mon. Not. Roy. Astron. Soc. , 2461(2013), arXiv:1309.3428 [astro-ph.HE].[26] J. Petrovi, P. D. Serpico, and G. Zaharijas, JCAP ,023 (2015), arXiv:1411.2980 [astro-ph.HE]. [27] Q. Yuan and B. Zhang, JHEAp , 1 (2014),arXiv:1404.2318 [astro-ph.HE].[28] Q. Yuan and K. Ioka, Astrophys. J. , 124 (2015),arXiv:1411.4363 [astro-ph.HE].[29] T. D. Brandt and B. Kocsis, Astrophys. J. , 15(2015), arXiv:1507.05616 [astro-ph.HE].[30] R. M. O’Leary, M. D. Kistler, M. Kerr, and J. Dexter,(2015), arXiv:1504.02477 [astro-ph.HE].[31] R. M. O’Leary, M. D. Kistler, M. Kerr, and J. Dexter,(2016), arXiv:1601.05797 [astro-ph.HE].[32] E. Carlson and S. Profumo, Phys.Rev. D90 , 023015(2014), arXiv:1405.7685 [astro-ph.HE].[33] D. Gaggero, M. Taoso, A. Urbano, M. Valli, and P. Ullio,JCAP , 056 (2015), arXiv:1507.06129 [astro-ph.HE].[34] I. Cholis, C. Evoli, F. Calore, T. Linden, C. Weniger, andD. Hooper, JCAP , 005 (2015), arXiv:1506.05119[astro-ph.HE].[35] E. Carlson, T. Linden, and S. Profumo, Phys. Rev.
D94 ,063504 (2016), arXiv:1603.06584 [astro-ph.HE].[36] M. Su, T. R. Slatyer, and D. P. Finkbeiner, Astrophys.J. , 1044 (2010), arXiv:1005.5480 [astro-ph.HE].[37] M. Ackermann et al. (Fermi-LAT), Astrophys. J. , 64(2014), arXiv:1407.7905 [astro-ph.HE].[38] J. Petrovi, P. D. Serpico, and G. Zaharija, JCAP ,052 (2014), arXiv:1405.7928 [astro-ph.HE].[39] M. Lisanti, S. Mishra-Sharma, N. L. Rodd, and B. R.Safdi, (2017), arXiv:1708.09385 [astro-ph.CO].[40] M. Lisanti, S. Mishra-Sharma, N. L. Rodd, B. R. Safdi,and R. H. Wechsler, (2017), arXiv:1709.00416 [astro-ph.CO].[41] A. Albert et al. (DES, Fermi-LAT), Astrophys. J. ,110 (2017), arXiv:1611.03184 [astro-ph.HE].[42] A. Geringer-Sameth, M. G. Walker, S. M. Koushiap-pas, S. E. Koposov, V. Belokurov, G. Torrealba, andN. W. Evans, Phys. Rev. Lett. , 081101 (2015),arXiv:1503.02320 [astro-ph.HE].[43] D. Hooper and T. Linden, JCAP , 016 (2015),arXiv:1503.06209 [astro-ph.HE].[44] A. Drlica-Wagner et al. (DES, Fermi-LAT), Astrophys.J. , L4 (2015), arXiv:1503.02632 [astro-ph.HE].[45] O. Macias, C. Gordon, R. M. Crocker, B. Coleman,D. Paterson, S. Horiuchi, and M. Pohl, (2016),arXiv:1611.06644 [astro-ph.HE].[46] R. Bartels, S. Krishnamurthy, and C. Weniger, Phys.Rev. Lett. , 051102 (2016), arXiv:1506.05104 [astro-ph.HE].[47] S. K. Lee, M. Lisanti, and B. R. Safdi, JCAP , 056(2015), arXiv:1412.6099 [astro-ph.CO].[48] S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer,and W. Xue, Phys. Rev. Lett. , 051103 (2016),arXiv:1506.05124 [astro-ph.HE].[49] T. Linden, N. L. Rodd, B. R. Safdi, and T. R. Slatyer,Phys. Rev.
D94 , 103013 (2016), arXiv:1604.01026 [astro-ph.HE].[50] S. Mishra-Sharma, N. L. Rodd, and B. R. Safdi, Astron.J. , 253 (2017), arXiv:1612.03173 [astro-ph.HE].[51] S. Horiuchi, M. Kaplinghat, and A. Kwa, JCAP ,053 (2016), arXiv:1604.01402 [astro-ph.HE].[52] F. Acero et al. (Fermi-LAT), Astrophys. J. Suppl. ,23 (2015), arXiv:1501.02003 [astro-ph.HE]. [53] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt,F. K. Hansen, M. Reinecke, and M. Bartelman, Astro-phys. J. , 759 (2005), arXiv:astro-ph/0409513 [astro-ph].[54] F. Yusef-Zadeh, J. Hewitt, M. Wardle, V. Tatis-cheff, D. Roberts, et al. , Astrophys.J. , 33 (2013),arXiv:1206.6882 [astro-ph.HE].[55] J. D. Simon et al. (DES), Astrophys. J. , 95 (2015),arXiv:1504.02889 [astro-ph.GA].[56] A. Geringer-Sameth, S. M. Koushiappas, and M. Walker,Astrophys. J. , 74 (2015), arXiv:1408.0002 [astro-ph.CO].[57] M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, et al. , JCAP , 051 (2011), arXiv:1012.4515 [hep-ph].[58] K. N. Abazajian and R. E. Keeley, Phys. Rev. D93 ,083514 (2016), arXiv:1510.06424 [hep-ph].[59] L. Zhang, H.-W. Rix, G. van de Ven, J. Bovy, C. Liu, andG. Zhao, Astrophys. J. , 108 (2013), arXiv:1209.0256[astro-ph.GA].[60] M. A. Snchez-Conde and F. Prada, Mon. Not. Roy.Astron. Soc. , 2271 (2014), arXiv:1312.1729 [astro-ph.CO]. [61] N. Klop, F. Zandanel, K. Hayashi, and S. Ando,Phys. Rev.
D95 , 123012 (2017), arXiv:1609.03509 [astro-ph.CO].[62] K. Hayashi, K. Ichikawa, S. Matsumoto, M. Ibe, M. N.Ishigaki, and H. Sugai, Mon. Not. Roy. Astron. Soc. ,2914 (2016), arXiv:1603.08046 [astro-ph.GA].[63] V. Bonnivard, C. Combet, D. Maurin, and M. G.Walker, Mon. Not. Roy. Astron. Soc. , 3002 (2015),arXiv:1407.7822 [astro-ph.HE].[64] Z. Ivezi´c, A. J. Connolly, J. T. VanderPlass, and A. Gray,
Statistics, Data Mining, and Machine Learning in As-tronomy: A Practical Python Guide for the Analysis ofSurvey Data (Princeton Series in Modern ObservationalAstronomy) (Princeton University Press, 2014).[65] M. Pato, F. Iocco, and G. Bertone, JCAP , 001(2015), arXiv:1504.06324 [astro-ph.GA].[66] C. F. McKee, A. Parravano, and D. J. Hollenbach, Astro-phys. J. , 13 (2015), arXiv:1509.05334 [astro-ph.GA].[67] A. W. McConnachie, Astron. J. , 4 (2012),arXiv:1204.1562 [astro-ph.CO].[68] M. Kaplinghat, S. Tulin, and H.-B. Yu, Phys. Rev. Lett.116