A deep search for decaying dark matter with XMM-Newton blank-sky observations
Joshua W. Foster, Marius Kongsore, Christopher Dessert, Yujin Park, Nicholas L. Rodd, Kyle Cranmer, Benjamin R. Safdi
AA deep search for decaying dark matter with
XMM-Newton blank-sky observations
Joshua W. Foster,
1, 2, 3, ∗ Marius Kongsore, Christopher Dessert,
1, 2, 3
YujinPark,
1, 2, 3
Nicholas L. Rodd,
2, 3
Kyle Cranmer, and Benjamin R. Safdi
2, 3, † Leinweber Center for Theoretical Physics, Department of Physics,University of Michigan, Ann Arbor, MI 48109, USA Berkeley Center for Theoretical Physics, University of California, Berkeley, CA 94720, USA Theoretical Physics Group, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA Center for Cosmology and Particle Physics, New York University, New York, NY 10003, USA
Sterile neutrinos with masses in the keV range are well-motivated extensions to the StandardModel that could explain the observed neutrino masses while also making up the dark matter (DM)of the Universe. If sterile neutrinos are DM then they may slowly decay into active neutrinosand photons, giving rise to the possibility of their detection through narrow spectral features inastrophysical X-ray data sets. In this work, we perform the most sensitive search to date for this andother decaying DM scenarios across the mass range from 5 to 16 keV using archival
XMM-Newton data. We reduce 547 Ms of data from both the MOS and PN instruments using observations takenacross the full sky and then use this data to search for evidence of DM decay in the ambient haloof the Milky Way. We determine the instrumental and astrophysical baselines with data taken faraway from the Galactic Center, and use Gaussian Process modeling to capture additional continuumbackground contributions. No evidence is found for unassociated X-ray lines, leading us to producethe strongest constraints to date on decaying DM in this mass range.
Sterile neutrino dark matter (DM) is a well-motivatedDM candidate that may give rise to observable nearlymonochromatic X-ray signatures [1–3]. In this scenariothe DM has a mass in the keV range and may decayinto an active neutrino and an X-ray, with energy setby half the rest mass of the sterile neutrino [4]. Sterileneutrino DM is motivated in part by the seesaw mecha-nism for explaining the active neutrino masses [5, 6]. Inthis work we present one of the most sensitive searchesfor sterile neutrino DM, along with other DM candidatesthat may decay to monochromatic X-rays, over the massrange m χ ∈ [5 ,
16] keV. We do so by searching for DMdecay from the ambient halo of the Milky Way using allarchival data from the
XMM-Newton telescope collectedfrom its launch until September 5, 2018.This work builds heavily off the method developed inDessert et al. [7], which used
XMM-Newton blank-sky ob-servations (BSOs) to strongly disfavor the decaying DMexplanation of the previously-observed 3.5 keV uniden-tified X-ray line (UXL). This UXL was found in nearbygalaxies and clusters [8–12]. However the analysis per-formed in Dessert et al. [7] was able to robustly rule outthe DM decay rate required to explain the previous 3.5keV UXL signals [13]. (For additional non-observations,see Refs. [14–20].) We extend the search in Dessert et al. [7] to the broader mass range m χ ∈ [5 , ∗ [email protected] † [email protected] Gaussian Process (GP) modeling to describe continuumresiduals, instead of parametric modeling as used in [7].As demonstrated in Dessert et al. [7], BSO searches forDM decaying in the Milky Way halo can be both moresensitive and more robust than extra-galactic searches,because (i) the expected DM flux, even at angles ∼ ◦ away from the GC, rivals the expected flux from the mostpromising extra-galactic objects, such as M31 and thePerseus cluster; (ii) promising extra-galactic targets havecontinuum and line-like X-ray features that are confound-ing backgrounds for DM searches (dwarf galaxies beingan exception [18, 21]), while BSOs instead focus on thelowest-background regions of the sky; (iii) extra-galactictargets require pointed observations, while in principleany observation collected by XMM-Newton is sensitiveto DM decay in the Milky Way, opening up considerablymore exposure time.The limits presented in this work represent thestrongest found using the
XMM-Newton instrument overthe energy range ∼ ∼ ◦ away in their case), while Ref. [26]searched for DM decay from M31 with NuSTAR.Our results put in tension efforts to explain the abun-dance of sterile neutrino DM in the Neutrino MinimalStandard Model ( ν MSM) [27–29], which may simultane-ously explain the observed neutrino masses, DM density,and baryon asymmetry. In this scenario the StandardModel is supplemented by three heavier sterile neutrinostates, the lightest of which is the DM candidate. TheDM abundance is generated through the mixing of sterileand active neutrinos [1], which can further be resonantlyenhanced by a finite lepton chemical potential [2, 29–35]. a r X i v : . [ a s t r o - ph . C O ] F e b r GC [deg] D - f a c t o r × − [ k e V / c m ] S i g n a l R O I Bkg. ROI D bkg Dark Matter Density Profile
MOS PN
Figure 1. Our fiducial D -factor, which is proportional to theexpected DM signal flux. Values are given in all 30 annuli,which are 6 ◦ wide in angular distance from the GC (with | b | > ◦ ), and we define a signal and background ROI asshown. In each ring, we compute the D -factor of all MOSor PN exposures, weighted according to the observation timeand field of view. The horizontal line indicates D bkg , themean D -factor in the background ROI. While sterile neutrino DM is one of the more motivatedscenarios that would produce a UXL, other models, suchas axion-like-particle DM [36], predict similar signatures.
Data reduction and processing.
We process and ana-lyze all publicly-available data collected before 5 Septem-ber 2018 by the metal oxide semiconductor (MOS) andpositive-negative (PN) cameras on board
XMM-Newton (see the Supplementary Material (SM) for a descriptionof our data-reduction pipeline). We subject each expo-sure to a set of quality cuts, which are described shortly.Those exposures satisfying the quality cuts are includedin our angularly-binned data products. In particular, wedivide the sky into 30 concentric annuli centered aroundthe GC, each with a width of 6 ◦ in angular radius fromthe GC, r GC , where cos( r GC ) = cos( l ) cos( b ) in termsof the Galactic longitude, l , and latitude, b . We labelthese from 1 to 30, starting from the innermost ring. Wefurther mask the Galactic Plane such that we only in-clude the region | b | ≥ ◦ . In each ring we then producestacked spectra where, in each energy bin, we sum overthe counts from each exposure whose central position lieswithin that annulus. We produce separate data sets forthe MOS and PN cameras, which have 2400 and 4096 en-ergy channels, respectively. In addition to stacking thecounts in each ring and energy channel, we also constructthe appropriately weighted detector response matrices inevery ring for forward modeling an incident astrophysicalflux. The full-sky maps and associated modeling data areprovided as Supplementary Data [37] in both the annuliand in finer-resolution HEALPix binning [38].We search for and mask point sources when process- ing the individual exposures, and for each exposure weconstruct the energy spectrum in the unmasked region.BSOs are expected to be dominated by the quiescentparticle background (QPB) caused by cosmic ray inter-actions with the detector. For regions of the sky thatare not focused on a bright Galactic or extra-galacticsource, the QPB counts should dominate over the extra-galactic X-ray background [39]. However, the QPB istime-dependent and will vary over exposures because of e.g. flaring activity (which, as described in the SM, wefilter for). We require that the 2 −
10 keV fluxes for theMOS (PN) exposures be less than the median across allobservations of 0.09 cts/keV/s (0.39 cts/keV/s) to be in-cluded in our analysis. Note that a putative signal wouldcontribute a negligible fraction to the 2 −
10 keV countrates and thus not change the selected observations [7].The cuts leave 215 Ms (57 Ms) of the total 438 Ms (109Ms) of full-sky ( | b | ≥ ◦ ) exposure time for MOS (PN).We analyze the MOS data from 2.5 to 8 keV and thePN data from 2.5 to 7 keV, in order to exclude intervalscontaining large instrumental features. Data analysis.
Having constructed our data in all 30rings, we divide the full sky into two regions of inter-est (ROI): a signal ROI, consisting of annuli 1 through8 (0 ◦ ≤ r GC ≤ ◦ with | b | ≥ ◦ ), inclusive, andthe background ROI, consisting of annuli 20 through 30(114 ◦ ≤ r GC ≤ ◦ with | b | ≥ ◦ ). The regions areillustrated in Fig. 1. The MOS (PN) exposure time inthe signal ROI is 25.27 Ms (5.56 Ms), whereas in thebackground ROI it is 62.51 Ms (17.54 Ms). The sig-nal flux is proportional to the D -factor, which is definedby the line-of-sight integral of the Galactic DM density ρ DM , D ≡ (cid:82) ds ρ DM . In Fig. 1 we show the appropriatelyweighted D -factor in each annuli. The motivation for thetwo ROIs is that the signal should dominate in the innerregions of the Galaxy and become progressively weakerfurther away from the GC. The background ROI is chosento be large enough to have significantly more exposuretime than the signal ROI, so that using the background-subtracted data does not significantly broaden the statis-tical uncertainties. We stack the data over the full back-ground ROI, which has D -factor D bkg , and use this as anestimate of the instrumental and astrophysical baselinefluxes by subtracting this data from the data in each ringof the signal ROI. This subtraction mostly removes largeinstrumental lines, as illustrated in SM Fig. S1.We analyze the background-subtracted data in eachannulus for evidence of a UXL. The data is modeledas a combination of narrow spectral features at the lo-cations of known astrophysical and instrumental lines,and a continuum flux which we account for using GPmodeling. Note that the instrumental lines need not becompletely removed by the data-subtraction procedure,leaving a residual flux or flux deficit that must be mod-eled. Astrophysical emission lines from the Milky Wayplasma should be brighter in the signal ROI, and so arealso expected to appear in the background-subtracteddata. For both astrophysical and instrumental lines, thelines are modeled using the forward modeling matricesfor MOS and PN. We allow the instrumental lines tohave either positive or negative normalizations, while theastrophysical lines are restricted to have positive nor-malizations. To decide which lines to include in ourresidual background model we start with an initial listof known instrumental and astrophysical lines. The in-strumental lines are determined from an analysis of thebackground ROI data, while the astrophysical lines arethose expected to be produced by the Milky Way. In eachring, and for MOS and PN independently, we then deter-mine the significance of each emission line, keeping thoseabove ∼ σ . As a result, every ring has a different set oflines included in the analysis. We note that it is conceiv-able that a UXL might be inadvertently removed by anoverly-subtracted instrumental line at the same energy;however, it would be highly unlikely for such a conspir-acy to occur in every ring, given the varying D -factor.The effects of sub-threshold instrumental lines are miti-gated through a spurious-signal nuisance parameter [40],as discussed in the SM.The unprecedented data volume incorporated intothis analysis necessitates a flexible approach to model-ing the residual continuum emission, which is accom-plished with GP modeling, in order to minimize back-ground mismodeling. As opposed to parametric mod-eling, where the model is specified by a specific func-tional form and associated list of model parameters,GP modeling is non-parametric: the model expecta-tions for the data at two different energies, E and E (cid:48) , are assumed to be normally distributed with non-trivial covariance. Taking the model expectation tohave zero mean, the GP model is then fully specifiedby the covariance kernel, K ( E, E (cid:48) ). We model themean-subtracted data using the non-stationary kernel K ( E, E (cid:48) ) = A GP exp (cid:2) − ( E − E (cid:48) ) / (2 EE (cid:48) σ E ) (cid:3) , imple-mented in george [41], where σ E is the correlation-lengthhyperparameter and A GP is the amplitude hyperparam-eter. We fix σ E such that it is larger than the energyresolution of the detector, which is δE/E ∼ .
03 acrossmost energies for MOS and PN, while ensuring σ E iskept small enough to have the flexibility to model realvariations in the data. The goal is to balance two com-peting effects. If σ E approaches the lower limit imposedby the energy resolution of the detector, then the GPmodel would have the flexibility to account for line-likefeatures, which would reduce our sensitivity when search-ing for such features over the continuum background. Onthe other hand, if σ E is too large then the GP continuummodel may not accurately model real small-scale varia-tions in the data. In our fiducial analysis we fix σ E = 0 . A GP is treated as a nuisance parameterthat is profiled over when searching for UXLs.We then follow the statistical approach developedin Frate et al. [42], which used GP modeling to perform . . . . . . Energy [keV] . . . . F l u x [ ph o t o n s / s / k e V ] t : 1.36 Null Fit Signal Fit
MOS Annulus 1
Figure 2. The background-subtracted MOS data for theinnermost annulus, downbinned by a factor of four for pre-sentation purposes. The indiciated best fit null and signalmodels, for a 3.5 keV UXL, are constructed using the GPmodeling described in the text. an improved search for narrow resonances over a con-tinuum background in the context of the Large HadronCollider. In particular, we construct a likelihood ratioΛ between the model with and without the signal com-ponent, where the signal is the UXL line at fixed energy E sig . The null model is as above, the combination of a GPmodel with a single nuisance parameter A GP , and a setof background lines, whose amplitudes are treated as nui-sance parameters. We use the marginal likelihood fromthe GP fit in the construction of the likelihood ratio [42].Note that as the number of counts in all energy binsis large ( (cid:29) t = − t follows a χ -distribution. The 95%one-sided upper limits are constructed from the profilelikelihood, as a function of the signal amplitude.We implement this procedure and scan for a UXL from2.5 to 8 keV in 5 eV intervals. At each test point we con-struct profile likelihoods for signal flux independently foreach ring using the background-subtracted MOS and PNdata. We then combine the likelihoods between rings– and eventually cameras – in a joint likelihood in thecontext of the DM model, as discussed shortly. As anexample, Fig. 2 illustrates the signal and null modelfits to the innermost MOS background-subtracted signal-annulus data for a putative UXL at 3.5 keV (indicated bythe vertical dashed line). Note that while the fit is per-formed over the full energy range (2.5 − σ and 2 σ uncertainties on the flux. In the bottom panel we il- − . − . . . . F l u x [ ph o t o n s / c m / s / s r ] MOS Annulus 13 4 5 6 7 8
Energy [keV] − S i g n i fi c a n ce ( √ t ) Figure 3. (Upper) The best-fit signal flux, and 1 and 2 σ uncertainties, as a function of the central UXL energy acrossour full energy range for the innermost MOS ring. (Lower)The corresponding significance in favor of the signal model,multiplied by the sign of the best fit UXL normalization atthat energy, along with the 1/2 σ expectations under the nullhypothesis. Included background lines are indicated in gray. lustrate the detection significance, √ t , multiplied by thesign of the best-fit flux, such that excesses and deficitsappear as positive and negative values, respectively. Thelocations of astrophysical and instrumental lines includedin the residual background model are indicated in grey,and we do not search for UXLs within 0.025 keV of these(as shaded). Note that we lose sensitivity in the vicinityof the background lines because of the degeneracy be-tween the signal and background-line spectral templates,leading to the uncertainties widening in Fig. 3 near theseenergies. DM interpretation.
We combine together the profilelikelihoods from the individual annuli to test the decay-ing DM model. In the context of sterile neutrino DMwith mass m χ and mixing angle θ , the DM decay in theGalactic halo produces an X-ray flux at energy m χ / ≈ .
26 photons / cm / s / sr × (cid:16) m χ . (cid:17) (cid:32) D keV / cm (cid:33) (cid:18) sin (2 θ )10 − (cid:19) , (1)Note that the D -factors, appropriately averaged overobservations in the individual annuli, are illustrated inFig. 1. Thus, at fixed DM mass m χ we may constructprofile likelihoods as functions of sin (2 θ ) by using (1)to appropriately combine the profile likelihoods as func-tions of flux in the individual annuli. We subtract off D bkg from the D -factors in each signal ring since anyUXL would also appear in the background ROI and thusbe included in the background subtraction. The D -factors may be computed from the DM densityprofile of the Milky Way. Modern hydrodynamic cosmo-logical simulations indicate that the DM density profilein Milky Way mass halos generally have a high degree ofspherical symmetry (for a review, see Ref. [44]). Further,the presence of baryons contracts the inner ∼
10 kpc ofthe profile away from the canonical Navarro, Frenk, andWhite (NFW) DM distribution [45, 46], so that there isan enhancement of the DM density at smaller radii ver-sus the NFW expectation [47–52], though cores could de-velop on top of this contraction at radii < ∼ Fire-2 simulations the DM-only and hydrodynamic sim-ulations produce DM density profiles that agree within ∼
25% at 10 kpc, but with baryons the density profilesare typically around twice as large as the NFW DM-onlyexpectation at distances ∼ ρ and a scale radius r s : ρ DM ( r ) = ρ / ( r/r s ) / (1+ r/r s ) .We use the recent results from Cautun et al. [57], whocombined Gaia
DR2 Galactic rotation curve data [58]with total mass estimates for the Galaxy from satelliteobservations [59, 60]. These data imply, in the context ofthe NFW model, a virial halo mass M DM200 = 0 . +0 . − . × M (cid:12) and a concentration c = r /r s = 13 . +3 . − . ,with a non-trivial covariance between M DM200 and c [57]such that lower concentrations prefer higher halo masses.Within the 2D 68% containment region for M DM200 and c quoted in Ref. [57], the lowest DM density at r ≈ . ρ = 6 . × M (cid:12) / kpc and r s =19 . ∼ .
29 GeV / cm , which isconsistent with local measurements of the DM densityusing the vertical motion of tracer stars perpendicular tothe Galactic plane, see, e.g. , Refs. [61, 62].We search for evidence of decaying DM in 10 eV inter-vals in mass between 5 −
16 keV, masking 0.1 keV win-dows around the locations of known lines, as indicated inFig. 4. We construct the joint likelihoods for the MOSand PN data sets. We test and account for additionalbackground mismodeling in the MOS and PN analysesby looking at the distribution of best-fit mixing angles inthe energy side-bands, using a technique similar to the“spurious signal” used by ATLAS in the search for theHiggs boson [40]. This procedure is described in the SMand only has a small effect at low masses. We then com-bine, at a given mass, the resulting MOS and PN profilelikelihoods to obtain the final profile likelihood used toconstruct the limit and discovery significance shown inFig. 4. In that figure we show the one-sided 95% up-per limit on sin (2 θ ) in the upper panel, along with the1 and 2 σ expectations for the power-constrained upperlimit [63] under the null hypothesis (shaded green and − − − − s i n ( θ ) DM under production (4)(3)(2)(1) m χ [keV] − S i g n i fi c a n ce ( √ t ) Figure 4. (Upper) The power-constrained 95% upper limiton the DM lifetime from this work, presented in the contextof the sterile-neutrino mixing angle sin (2 θ ), as a functionof the DM mass m χ . The dark grey regions correspond totheoretical bounds for DM underproduction in the ν MSM orbounds from previous X-ray searches (1)–(5); see text for de-tails. (Lower) The associated sign-weighted significance forthe UXL. Vertical grey regions denote background lines andare at least partially masked. Green and gold regions in-dicate 1/2 σ expectations under the null hypothesis. Theseresults are shown in the context of more general DM modelsas constraints on the DM lifetime in SM Fig. S6. gold, respectively).We find no evidence for decaying DM signals above ourpre-determined significance threshold of 5 σ global signifi-cance (corresponding to ∼ σ local significance), as shownin the bottom panel. In that figure we compare our upperlimit to previous limits in the literature, adjusted to ourfiducial DM model for the Milky Way where appropriate.In the context of the ν MSM it is impossible to explainall of the observed DM in the region marked “DM un-der production” because of the big bang nucleosynthesisbound on the lepton chemical potential [30–32]. Notethat the ν MSM also predicts that the DM becomes in-creasingly warm for decreasing m χ , which leads to ten-sion with Milky Way satellite galaxy counts for low m χ :data from the Dark Energy Survey and other Galacticsatellite surveys [64] constrains m χ greater than ∼ ν MSM [65] (which can be strengthened furtherwhen combined with strong lensing measurements [66]),though we note that our results apply to more generalDM production mechanisms that do not predict modifi-cations to small-scale structure. In Fig. 4 we also showprevious X-ray limits from (1) Dessert et al. [7], (2) a
Chandra search for DM decay in the Milky Way [67],(3) a
Chandra search for DM decay in M31 [14], and (4)combined
NuSTAR searches for DM decay: in the MilkyWay [22–24], the Bullet Cluster [25], and M31 [26]. Note that the results from Milky Way searches have been ad-justed to use the same DM density profile as in our fidu-cial analysis. In the SM we perform multiple systematictests of our analysis framework, such as showing that syn-thetic signals are correctly recovered, that point to theresults in Fig. 4 robustly ruling out mixing angles abovethe indicated curve.
Discussion.
In this work we search for evidence of de-caying DM to narrow X-ray lines for masses m χ ∈ [5 , XMM-Newton data set, accumulated from nearly 20 years ofobservations. The overarching strategy behind this workis that all observations taken by the telescope necessar-ily stare through a column of the Milky Way’s DM halo;by masking point sources in the images, we are able toreduce 547 Ms – more than 17 years – of data across thefull sky that can be used for a search for decaying DMin the Milky Way’s halo. However, incorporating thisunprecedented volume of data into a statistical analy-sis necessitates a careful consideration of the continuummodels and possible systematic uncertainties. We adopta data-driven approach, whereby we use data from faraway from the GC as a measure of the null hypothesisand also incorporate GP modeling to account for residualcontinuum emission. We find no significant evidence fordecaying DM, which leads us to set some of the strongestconstraints on the DM lifetime.Given the data volume incorporated into this analysisit is unlikely that further analyses of
XMM-Newton data,or
Chandra data, could produce qualitatively stronger re-sults on the DM lifetime in the mass range consideredhere. However, the approach taken in this work maylead to a powerful advancement in discovery power withfuture data sets from surveys such as those by the up-coming
Athena [68] and
XRISM [69] telescopes. A com-bination of the data collected by those missions and theanalysis framework introduced in this work may lead tothe discovery of decaying DM in the few-keV mass rangeat lifetimes beyond those probed in this work.
Acknowledgments.
We thank Kerstin Perez andChristoph Weniger for useful conversations. This workwas supported in part by the DOE Early Career GrantDESC0019225. This research used resources of theNational Energy Research Scientific Computing Cen-ter (NERSC) and the Lawrencium computational clus-ter provided by the IT Division at the Lawrence Berke-ley National Laboratory, supported by the Director, Of-fice of Science, and Office of Basic Energy Sciences, ofthe U.S. Department of Energy under Contract No. DE-AC02-05CH11231. NLR is supported by the Miller In-stitute for Basic Research in Science at the University ofCalifornia, Berkeley. KC is partially supported by NSFgrant PHY-1505463m, NSF awards ACI-1450310, OAC-1836650, and OAC-1841471, and the Moore-Sloan DataScience Environment at NYU. [1] Scott Dodelson and Lawrence M. Widrow, “Sterile-neutrinos as dark matter,” Phys. Rev. Lett. , 17–20(1994), arXiv:hep-ph/9303287 [hep-ph].[2] Xiang-Dong Shi and George M. Fuller, “A New dark mat-ter candidate: Nonthermal sterile neutrinos,” Phys. Rev.Lett. , 2832–2835 (1999), arXiv:astro-ph/9810076.[3] Alexander Kusenko, “Sterile neutrinos, dark matter,and the pulsar velocities in models with a Higgs sin-glet,” Phys. Rev. Lett. , 241301 (2006), arXiv:hep-ph/0609081 [hep-ph].[4] Palash B. Pal and Lincoln Wolfenstein, “Radiative decaysof massive neutrinos,” Phys. Rev. D , 766–773 (1982).[5] Tsutomu Yanagida, “Horizontal Symmetry and Massesof Neutrinos,” Prog. Theor. Phys. , 1103 (1980).[6] Rabindra N. Mohapatra and Goran Senjanovic, “Neu-trino Mass and Spontaneous Parity Violation,” Phys.Rev. Lett. , 912 (1980), [,231(1979)].[7] Christopher Dessert, Nicholas L. Rodd, and Benjamin R.Safdi, “The dark matter interpretation of the 3.5-keVline is inconsistent with blank-sky observations,” Science , 1465 (2020), arXiv:1812.06976 [astro-ph.CO].[8] Esra Bulbul, Maxim Markevitch, Adam Foster, Ran-dall K. Smith, Michael Loewenstein, and Scott W. Ran-dall, “Detection of An Unidentified Emission Line in theStacked X-ray spectrum of Galaxy Clusters,” Astrophys.J. , 13 (2014), arXiv:1402.2301 [astro-ph.CO].[9] Alexey Boyarsky, Oleg Ruchayskiy, Dmytro Iakubovskyi,and Jeroen Franse, “Unidentified Line in X-Ray Spectraof the Andromeda Galaxy and Perseus Galaxy Cluster,”Phys. Rev. Lett. , 251301 (2014), arXiv:1402.4119[astro-ph.CO].[10] O. Urban, N. Werner, S. W. Allen, A. Simionescu, J. S.Kaastra, and L. E. Strigari, “A Suzaku Search for DarkMatter Emission Lines in the X-ray Brightest GalaxyClusters,” Mon. Not. Roy. Astron. Soc. , 2447–2461(2015), arXiv:1411.0050 [astro-ph.CO].[11] Tesla E. Jeltema and Stefano Profumo, “Discovery of a3.5 keV line in the Galactic Centre and a critical lookat the origin of the line across astronomical targets,”Mon. Not. Roy. Astron. Soc. , 2143–2152 (2015),arXiv:1408.1699 [astro-ph.HE].[12] Nico Cappelluti, Esra Bulbul, Adam Foster, PriyamvadaNatarajan, Megan C. Urry, Mark W. Bautz, FrancescaCivano, Eric Miller, and Randall K. Smith, “Searchingfor the 3.5 keV Line in the Deep Fields with Chandra:the 10 Ms observations,” Astrophys. J. , 179 (2018),arXiv:1701.07932 [astro-ph.CO].[13] Christopher Dessert, Nicholas L. Rodd, and Ben-jamin R. Safdi, “Response to a comment on Dessert etal. ”The dark matter interpretation of the 3.5 keV lineis inconsistent with blank-sky observations”,” (2020),arXiv:2006.03974 [astro-ph.CO].[14] Shunsaku Horiuchi, Philip J. Humphrey, Jose Onorbe,Kevork N. Abazajian, Manoj Kaplinghat, and SheaGarrison-Kimmel, “Sterile neutrino dark matter boundsfrom galaxies of the Local Group,” Phys. Rev. D89 ,025017 (2014), arXiv:1311.0282 [astro-ph.CO].[15] D. Malyshev, A. Neronov, and D. Eckert, “Constraintson 3.55 keV line emission from stacked observationsof dwarf spheroidal galaxies,” Phys. Rev.
D90 , 103506(2014), arXiv:1408.3531 [astro-ph.HE]. [16] Michael E. Anderson, Eugene Churazov, and Joel N.Bregman, “Non-Detection of X-Ray Emission From Ster-ile Neutrinos in Stacked Galaxy Spectra,” Mon. Not. Roy.Astron. Soc. , 3905–3923 (2015), arXiv:1408.4115[astro-ph.HE].[17] Takayuki Tamura, Ryo Iizuka, Yoshitomo Maeda,Kazuhisa Mitsuda, and Noriko Y. Yamasaki, “An X-ray Spectroscopic Search for Dark Matter in the PerseusCluster with Suzaku,” Publ. Astron. Soc. Jap. , 23(2015), arXiv:1412.1869 [astro-ph.HE].[18] Tesla E. Jeltema and Stefano Profumo, “Deep XMMObservations of Draco rule out at the 99% ConfidenceLevel a Dark Matter Decay Origin for the 3.5 keV Line,”Mon. Not. Roy. Astron. Soc. , 3592–3596 (2016),arXiv:1512.01239 [astro-ph.HE].[19] F. A. Aharonian et al. (Hitomi), “ Hitomi constraintson the 3.5 keV line in the Perseus galaxy cluster,” As-trophys. J. , L15 (2017), arXiv:1607.07420 [astro-ph.HE].[20] A. Gewering-Peine, D. Horns, and J. H. M. M. Schmitt,“A sensitive search for unknown spectral emission linesin the diffuse X-ray background with XMM-Newton,”JCAP , 036 (2017), arXiv:1611.01733 [astro-ph.HE].[21] Oleg Ruchayskiy, Alexey Boyarsky, Dmytro Iakubovskyi,Esra Bulbul, Dominique Eckert, Jeroen Franse, DenysMalyshev, Maxim Markevitch, and Andrii Neronov,“Searching for decaying dark matter in deep XMM–Newton observation of the Draco dwarf spheroidal,”Mon. Not. Roy. Astron. Soc. , 1390–1398 (2016),arXiv:1512.07217 [astro-ph.HE].[22] A. Neronov, Denys Malyshev, and Dominique Eck-ert, “Decaying dark matter search with NuSTAR deepsky observations,” Phys. Rev. D , 123504 (2016),arXiv:1607.07328 [astro-ph.HE].[23] Kerstin Perez, Kenny C. Y. Ng, John F. Beacom, CoraHersh, Shunsaku Horiuchi, and Roman Krivonos, “Al-most closing the ν MSM sterile neutrino dark matter win-dow with NuSTAR,” Phys. Rev. D , 123002 (2017),arXiv:1609.00667 [astro-ph.HE].[24] Brandon M. Roach, Kenny C.Y. Ng, Kerstin Perez,John F. Beacom, Shunsaku Horiuchi, Roman Krivonos,and Daniel R. Wik, “NuSTAR Tests of Sterile-NeutrinoDark Matter: New Galactic Bulge Observations andCombined Impact,” Phys. Rev. D , 103011 (2020),arXiv:1908.09037 [astro-ph.HE].[25] S. Riemer-Sørensen et al. , “Dark matter line emissionconstraints from NuSTAR observations of the BulletCluster,” Astrophys. J. , 48 (2015), arXiv:1507.01378[astro-ph.CO].[26] Kenny C.Y. Ng, Brandon M. Roach, Kerstin Perez,John F. Beacom, Shunsaku Horiuchi, Roman Krivonos,and Daniel R. Wik, “New Constraints on Sterile Neu-trino Dark Matter from NuST AR
M31 Observations,”Phys. Rev. D , 083005 (2019), arXiv:1901.01262 [astro-ph.HE].[27] Takehiko Asaka, Steve Blanchet, and Mikhail Sha-poshnikov, “The nuMSM, dark matter and neutrinomasses,” Phys. Lett. B , 151–156 (2005), arXiv:hep-ph/0503065.[28] Takehiko Asaka and Mikhail Shaposhnikov, “The ν MSM,dark matter and baryon asymmetry of the universe,”
Phys. Lett. B , 17–26 (2005), arXiv:hep-ph/0505013.[29] Laurent Canetti, Marco Drewes, Tibor Frossard, andMikhail Shaposhnikov, “Dark Matter, Baryogenesis andNeutrino Oscillations from Right Handed Neutrinos,”Phys. Rev. D , 093006 (2013), arXiv:1208.4607 [hep-ph].[30] A.D. Dolgov, S.H. Hansen, S. Pastor, S.T. Petcov,G.G. Raffelt, and D.V. Semikoz, “Cosmological boundson neutrino degeneracy improved by flavor oscilla-tions,” Nucl. Phys. B , 363–382 (2002), arXiv:hep-ph/0201287.[31] Pasquale D. Serpico and Georg G. Raffelt, “Lepton asym-metry and primordial nucleosynthesis in the era of pre-cision cosmology,” Phys. Rev. D , 127301 (2005),arXiv:astro-ph/0506162.[32] Alexey Boyarsky, Oleg Ruchayskiy, and Mikhail Sha-poshnikov, “The Role of sterile neutrinos in cosmologyand astrophysics,” Ann. Rev. Nucl. Part. Sci. , 191–214 (2009), arXiv:0901.0011 [hep-ph].[33] M. Laine and M. Shaposhnikov, “Sterile neutrino darkmatter as a consequence of nuMSM-induced lepton asym-metry,” JCAP , 031 (2008), arXiv:0804.4543 [hep-ph].[34] Tejaswi Venumadhav, Francis-Yan Cyr-Racine,Kevork N. Abazajian, and Christopher M. Hirata,“Sterile neutrino dark matter: Weak interactions in thestrong coupling epoch,” Phys. Rev. D , 043515 (2016),arXiv:1507.06655 [astro-ph.CO].[35] John F. Cherry and Shunsaku Horiuchi, “Closing in onResonantly Produced Sterile Neutrino Dark Matter,”Phys. Rev. D , 083015 (2017), arXiv:1701.07874 [hep-ph].[36] Tetsutaro Higaki, Kwang Sik Jeong, and FuminobuTakahashi, “The 7 keV axion dark matter and theX-ray line signal,” Phys. Lett. B , 25–31 (2014),arXiv:1402.6965 [hep-ph].[37] See /github.com/bsafdi/XMM BSO DATA for the list ofobservations used and processed data products.[38] K.M. Gorski, Eric Hivon, A.J. Banday, B.D. Wan-delt, F.K. Hansen, M. Reinecke, and M. Bartelman,“HEALPix - A Framework for high resolution dis-cretization, and fast analysis of data distributed on thesphere,” Astrophys. J. , 759–771 (2005), arXiv:astro-ph/0409513.[39] D. H. Lumb, R. S. Warwick, M. Page, and A. De Luca,“X-ray background measurements with xmm-newtonepic,” Astron. Astrophys. , 93 (2002), arXiv:astro-ph/0204147 [astro-ph].[40] Georges Aad et al. (ATLAS), “Measurement of Higgs bo-son production in the diphoton decay channel in pp col-lisions at center-of-mass energies of 7 and 8 TeV withthe ATLAS detector,” Phys. Rev. D , 112015 (2014),arXiv:1408.7084 [hep-ex].[41] S. Ambikasaran, D. Foreman-Mackey, L. Greengard,D. W. Hogg, and M. O’Neil, “Fast Direct Methods forGaussian Processes,” (2014).[42] Meghan Frate, Kyle Cranmer, Saarik Kalia, Alexan-der Vandenberg-Rodes, and Daniel Whiteson, “Mod-eling Smooth Backgrounds and Generic Localized Sig-nals with Gaussian Processes,” (2017), arXiv:1709.05681[physics.data-an].[43] Palash B. Pal and Lincoln Wolfenstein, “Radiative De-cays of Massive Neutrinos,” Phys. Rev. D , 766 (1982).[44] Mark Vogelsberger, Federico Marinacci, Paul Torrey,and Ewald Puchwein, “Cosmological Simulations of Galaxy Formation,” Nature Rev. Phys. , 42–66 (2020),arXiv:1909.07976 [astro-ph.GA].[45] Julio F. Navarro, Carlos S. Frenk, and Simon D. M.White, “The Structure of cold dark matter halos,” As-trophys. J. , 563–575 (1996), astro-ph/9508025.[46] Julio F. Navarro, Carlos S. Frenk, and Simon D. M.White, “A Universal density profile from hierarchi-cal clustering,” Astrophys. J. , 493–508 (1997),arXiv:astro-ph/9611107 [astro-ph].[47] Oleg Y. Gnedin, Andrey V. Kravtsov, Anatoly A. Klypin,and Daisuke Nagai, “Response of dark matter halos tocondensation of baryons: Cosmological simulations andimproved adiabatic contraction model,” Astrophys. J. , 16–26 (2004), arXiv:astro-ph/0406247.[48] Matthieu Schaller, Carlos S. Frenk, Richard G. Bower,Tom Theuns, Adrian Jenkins, Joop Schaye, Robert A.Crain, Michelle Furlong, Claudio Dalla Vecchia, andI.G. McCarthy, “Baryon effects on the internal struc-ture of ΛCDM haloes in the EAGLE simulations,”Mon. Not. Roy. Astron. Soc. , 1247–1267 (2015),arXiv:1409.8617 [astro-ph.CO].[49] Qirong Zhu, Federico Marinacci, Moupiya Maji, YuexingLi, Volker Springel, and Lars Hernquist, “Baryonic im-pact on the dark matter distribution in Milky Way-sizedgalaxies and their satellites,” Mon. Not. Roy. Astron. Soc. , 1559–1580 (2016), arXiv:1506.05537 [astro-ph.CO].[50] Aaron A. Dutton, Andrea V. Macci`o, Avishai Dekel,Liang Wang, Gregory S. Stinson, Aura Obreja, AriannaDi Cintio, Chris B. Brook, Tobias Buck, and Xi Kang,“NIHAO IX: The role of gas inflows and outflows in driv-ing the contraction and expansion of cold dark matterhaloes,” Mon. Not. Roy. Astron. Soc. , 2658–2675(2016), arXiv:1605.05323 [astro-ph.GA].[51] Philip F Hopkins et al. , “FIRE-2 Simulations: Physicsversus Numerics in Galaxy Formation,” Mon. Not. Roy.Astron. Soc. , 800–863 (2018), arXiv:1702.06148[astro-ph.GA].[52] Mark R. Lovell et al. , “The fraction of dark mat-ter within galaxies from the IllustrisTNG simulations,”Mon. Not. Roy. Astron. Soc. , 1950–1975 (2018),arXiv:1801.10170 [astro-ph.GA].[53] T.K. Chan, D. Kereˇs, J. O˜norbe, P.F. Hopkins, A.L. Mu-ratov, C. A. Faucher-Gigu`ere, and E. Quataert, “Theimpact of baryonic physics on the structure of dark mat-ter haloes: the view from the FIRE cosmological sim-ulations,” Mon. Not. Roy. Astron. Soc. , 2981–3001(2015), arXiv:1507.02282 [astro-ph.GA].[54] Pol Mollitor, Emmanuel Nezri, and Romain Teyssier,“Baryonic and dark matter distribution in cosmologi-cal simulations of spiral galaxies,” Mon. Not. Roy. As-tron. Soc. , 1353–1369 (2015), arXiv:1405.4318 [astro-ph.GA].[55] Matthieu Portail, Ortwin Gerhard, Christopher Wegg,and Melissa Ness, “Dynamical modelling of the galac-tic bulge and bar: the Milky Way’s pattern speed, stel-lar and dark matter mass distribution,” Mon. Not. Roy.Astron. Soc. , 1621–1644 (2017), arXiv:1608.07954[astro-ph.GA].[56] Alexandres Lazar, James S. Bullock, Michael Boylan-Kolchin, T. K. Chan, Philip F. Hopkins, Andrew S.Graus, Andrew Wetzel, Kareem El-Badry, CoralWheeler, Maria C. Straight, Duˇsan Kereˇs, Claude-Andr´e Faucher-Gigu`ere, Alex Fitts, and Shea Garrison-Kimmel, “A dark matter profile to model diverse feedback-induced core sizes of ΛCDM haloes,” MNRAS , 2393–2417 (2020), arXiv:2004.10817 [astro-ph.GA].[57] Marius Cautun, Alejandro Ben´ıtez-Llambay, Alis J.Deason, Carlos S. Frenk, Azadeh Fattahi, Facundo A.G´omez, Robert J. J. Grand, Kyle A. Oman, Julio F.Navarro, and Christine M. Simpson, “The milky waytotal mass profile as inferred from Gaia DR2,” MNRAS , 4291–4313 (2020), arXiv:1911.04557 [astro-ph.GA].[58] Anna-Christina Eilers, David W. Hogg, Hans-WalterRix, and Melissa K. Ness, “The Circular Velocity Curveof the Milky Way from 5 to 25 kpc,” Astrophys. J. ,120 (2019), arXiv:1810.09466 [astro-ph.GA].[59] Laura L. Watkins, Roeland P. van der Marel,Sangmo Tony Sohn, and N. Wyn Evans, “Evidencefor an Intermediate-mass Milky Way from Gaia DR2Halo Globular Cluster Motions,” ApJ , 118 (2019),arXiv:1804.11348 [astro-ph.GA].[60] Thomas M. Callingham, Marius Cautun, Alis J. Dea-son, Carlos S. Frenk, Wenting Wang, Facundo A.G´omez, Robert J. J. Grand, Federico Marinacci, andRuediger Pakmor, “The mass of the Milky Way fromsatellite dynamics,” MNRAS , 5453–5467 (2019),arXiv:1808.10456 [astro-ph.GA].[61] J.I. Read, “The Local Dark Matter Density,” J. Phys. G , 063101 (2014), arXiv:1404.1938 [astro-ph.GA].[62] Pablo F. de Salas and Axel Widmark, “Dark matter lo-cal density determination: recent observations and futureprospects,” (2020), arXiv:2012.11477 [astro-ph.GA].[63] Glen Cowan, Kyle Cranmer, Eilam Gross, andOfer Vitells, “Power-Constrained Limits,” (2011),arXiv:1105.3166 [physics.data-an].[64] A. Drlica-Wagner et al. (DES), “Milky Way SatelliteCensus. I. The Observational Selection Function forMilky Way Satellites in DES Y3 and Pan-STARRSDR1,” Astrophys. J. , 1 (2020), arXiv:1912.03302[astro-ph.GA].[65] E. O. Nadler et al. (DES), “Milky Way Satellite Cen-sus. III. Constraints on Dark Matter Properties fromObservations of Milky Way Satellite Galaxies,” (2020),arXiv:2008.00022 [astro-ph.CO].[66] Ethan O. Nadler, Simon Birrer, Daniel Gilman, Risa H.Wechsler, Xiaolong Du, Andrew Benson, Anna M.Nierenberg, and Tommaso Treu, “Dark Matter Con-straints from a Unified Analysis of Strong GravitationalLenses and Milky Way Satellite Galaxies,” (2021),arXiv:2101.07810 [astro-ph.CO].[67] Dominic Sicilian, Nico Cappelluti, Esra Bulbul,Francesca Civano, Massimo Moscetti, and Christo-pher S. Reynolds, “Probing the Milky Way’s Dark Mat-ter Halo for the 3.5 keV Line,” Astrophys. J. , 146(2020), arXiv:2008.02283 [astro-ph.HE].[68] X. Barcons, K. Nandra, D. Barret, J.W. den Herder,A.C. Fabian, L. Piro, and M.G. Watson (Athena Team),“Athena: the X-ray observatory to study the hot and en-ergetic Universe,” J. Phys. Conf. Ser. , 012008 (2015).[69] “Science with the X-ray Imaging and Spectroscopy Mis-sion (XRISM),” (2020), arXiv:2003.04962 [astro-ph.HE].[70] “Users Guide to the XMM-Newton Science Analysis Sys-tem”, Issue 14.0, 2018 (ESA: XMM-Newton SOC).[71] Glen Cowan, Kyle Cranmer, Eilam Gross, and OferVitells, “Asymptotic formulae for likelihood-based testsof new physics,” Eur. Phys. J. C71 , 1554 (2011), [Er-ratum: Eur. Phys. J.C73,2501(2013)], arXiv:1007.1727[physics.data-an]. [72] de Plaa, J., Werner, N., Simionescu, A., Kaastra, J. S.,Grange, Y. G., and Vink, J., “Cold fronts and multi-temperature structures in the core of abell 2052,” A&A , A81 (2010).[73] Leccardi, A. and Molendi, S., “Radial temperature pro-files for a large sample of galaxy clusters observed withxmm-newton*,” A&A , 359–373 (2008).[74] Adam Foster, Randall K. Smith, Nancy S. Brickhouse,and Xiaohong Cui, “AtomDB and PyAtomDB: AtomicData and Modelling Tools for High Energy and Non-Maxwellian Plasmas,” in
American Astronomical SocietyMeeting Abstracts , American Astronomical SocietyMeeting Abstracts, Vol. 227 (2016) p. 211.08.[75] Katsuji Koyama et al. , “Iron and Nickel Line Diagnosticsfor the Galactic Center Diffuse Emission,” Publ. Astron.Soc. Jap. , 245 (2007), arXiv:astro-ph/0609215.[76] M. Ackermann et al. (Fermi-LAT), “Search for Gamma-ray Spectral Lines with the Fermi Large Area Telescopeand Dark Matter Implications,” Phys. Rev. D , 082002(2013), arXiv:1305.5597 [astro-ph.HE].[77] Andrea Albert, German A. Gomez-Vargas, MichaelGrefe, Carlos Munoz, Christoph Weniger, Elliott D.Bloom, Eric Charles, Mario N. Mazziotta, and AldoMorselli (Fermi-LAT), “Search for 100 MeV to 10 GeV γ -ray lines in the Fermi-LAT data and implications forgravitino dark matter in µν SSM,” JCAP , 023 (2014),arXiv:1406.3430 [astro-ph.HE].[78] M. Ackermann et al. (Fermi-LAT), “Updated search forspectral lines from Galactic dark matter interactionswith pass 8 data from the Fermi Large Area Telescope,”Phys. Rev. D , 122002 (2015), arXiv:1506.00013 [astro-ph.HE]. Supplementary Material for A deep search for decaying dark matter with
XMM-Newton blank-sky observations
Joshua W. Foster, Marius Kongsore, Christopher Dessert, Yujin Park, Nicholas L. Rodd, Kyle Cranmer, andBenjamin R. SafdiThis Supplementary Material (SM) is organized as follows. In Sec. I we provide additional details behind our datareduction and analysis procedure. In Sec. II we provide additional results from the main analysis presented in thisLetter. Sec. III presents non-trivial checks of our analysis procedure using synthetic signals. Lastly, in Sec. IV weperform multiple analysis variations to demonstrate the robustness of our main results.
I. DATA REDUCTION AND ANALYSIS
In this section, we detail our process for data reduction and analysis.
A. Data Reduction
We selected all XMM-Newton observations performed until September 5, 2018. For each of these observations, weretrieved the raw data products from the XMM-Newton Science Archive. For data reduction, we used the XMM-Newton Extended Source Analysis Software (ESAS) package, which is a part of the Science Analysis System [70] (SAS)version 17.0, and used for modeling sources covering the full XMM-Newton field-of-view and diffuse backgrounds.The data reduction process is described in detail in Ref. [7]; here, we summarize the important steps and pointout any differences. To reduce a given observation, we obtain the list of science exposures (i.e. pointings taken in amode usable for scientific purposes) from the summary files. For each exposure (independent of camera), we generatean event list and filter this list to only include events which were recorded during a period of low-background, whichcuts contamination from soft-proton flares. We then mask point sources in the field of view which contribute in anyenergy range (c.f. Ref. [7] where we only masked point sources in the 3-4 keV range). We also mask data from CCDsoperating in anomalous states. From the filtered and masked data products we create the photon-count data, theancillary response file (ARF), and the redistribution matrix file (RMF).The reduced data contains 11,805 observations, with 21,388 and 8,190 individual MOS and PN exposures, totalling438 Ms and 109 Ms of data. Given our focus is on searching for DM emission in otherwise dark regions of the sky,we place a cut on these data sets to isolate the astrophysically quietest amongst them. In particular, we constructthe integrated flux from 2 −
10 keV in all exposures, and determine the median value for MOS and PN separately as0.09 photons/keV/s and 0.39 photons/keV/s. All observations with integrated fluxes higher than these median valuesare excluded. This cut will remove observations with above average astrophysical emission, but also those wherethere is large instrumental or QPB emission (c.f. Ref. [7] where a separate cut on the QPB emission was performed).Further, we emphasize that even in the most optimistic scenario, a DM UXL will only provide an exceptionally smallcontribution to the total integrated flux, and thus this cut will not bias against a potential signal. In addition toremoving these bright exposures, we place two additional cuts. Firstly, all exposures with less than 500s of data areremoved, as the flux in such short exposures can be poorly characterized. Finally, we exclude all observations within2 ◦ of the plane of the Milky Way, which excises only a small amount of the expected DM signal, but a much largerfraction of the expected astrophysical emission associated with emission from our own galaxy.Exposure passing all three cuts are then divided into 30 rings, each of width 6 ◦ from the GC as described in the mainbody. The rings, numbered 1-30 starting from the GC, are used to form our signal ROI (rings 1-8) and backgroundROI (rings 20-30). B. Public Data Products
The processed data used to perform the analysis in this work is made fully publicly available atgithub.com/bsafdi/XMM BSO DATA. There we provide all the data required to reproduce our results. In par-ticular, we provide the data after the cuts described above in each ring for the MOS and PN cameras separately. Theinstrument response files, appropriately weighted across the exposures in each ring, are provided.
C. Analysis
In this section we provide additional details behind the analysis framework used to interpret the data productsdescribed above in the context of the decaying DM model. First, we describe how we analyze the flux data in theindividual annuli, and then we detail how those results are joined together to constrain the DM lifetime. Lastly, wedescribe how we test for and incorporate systematic uncertainties.
1. Construction of the profile likelihood
Let us first focus on the analysis of the (either MOS or PN) data in an individual ring k ∈ [1 , d k in this ring consists of background subtracted count rates d ki in each energy channel i . The count rates haveunits of cts / s / keV, as illustrated in e.g. Fig. S1, with Poisson counting uncertainties σ ki that arise from combiningthe statistical uncertainties in the signal and background data sets in the large-count limit, where the uncertaintiesbecome normally distributed. Our goal is then to compute the log-likelihood log p ( d k | θ ) as a function of the modelparameters θ = { A sig , θ nuis } , which consist of our parameter of interest, A sig , and our nuisance parameters θ nuis .The nuisance parameters include background line amplitudes, A j , with j indexing the different lines at energies E j ,and also hyperparameters for the GP model. In our fiducial analysis the only GP model hyperparameter is theamplitude of the double-exponential kernel A GP . Note that, as described shortly, in determining the instrumentalline list we also assign nuisance parameters to the locations of the lines. Our goal is to construct the profile likelihoodlog p ( d k | A sig ) = max nuis log p ( d k | θ ).Before describing the log-likelihood function in detail, we note that because we are in the large-count limit, so thatthe statistical fluctuations are normally distributed, we may interchangeably use the concept of modeling the data asthe sum of model components and subtracting model components from the data and considering the residuals. In thesmall-count limit, where the Poisson fluctuations are not nearly Gaussian, this approach would not be appropriate.The log-likelihood function that we use is a modification of the zero-mean GP marginal likelihood. The modificationthat we implement incorporates the background lines and the signal line of interest. For a given set of model parameters { A sig , A j } we construct the modified data vector y ki ( θ ) ≡ d ki − A sig µ k sig ,i − (cid:88) j A j µ kj,i − (cid:10) d ki − A sig µ k sig ,i − (cid:88) j A j µ kj,i (cid:11) i , (S1)where µ sig is the spectral template of the signal line of interest, with fixed normalization, as obtained by appropriatelysumming the forward modeling matrices of the individual exposures that compose the observations within the ringof interest, k . Similarly, µ kj,i denotes the fixed-normalization spectral template of the j th background line, at energy E j , in ring k (recall that i labels the detector energy channel). The quantity (cid:104)· · · (cid:105) i in (S1) denotes the average overenergy bins i , which implies that by construction the y ki ( θ ) have zero mean when averaged over the full energy rangeof the analysis. We postulate that the y ki are described by GP models, so that we may use the GP marginal likelihoodto compute the hyperparameter A GP :log p ( d k | θ ) = − y kT (cid:2) K + ( σ k ) I (cid:3) − y k −
12 log | K + ( σ k ) I | − n π ) . (S2)Above, n is the number of energy channels, and all matrix operations are taken in the space of energy channels, with( σ k ) I denoting the diagonal matrix with entries ( σ ki ) . The matrix K denotes the GP kernel. We implement thenon-stationary kernel K ( E, E (cid:48) ) = A GP exp (cid:20) − ( E − E (cid:48) ) EE (cid:48) σ E (cid:21) , (S3)which has the hyperparameters A GP and σ E . Note that later in the SM we show that similar results are obtainedusing the more standard double exponential kernel, but we chose the form of the kernel in (S3) for reasons discussedbelow.It is worth emphasizing that we have made the choice to describe the residuals of the background-subtracted data,after also subtracting the contributions from the instrumental lines, by a zero-mean GP model. An alternative strategy Note that the line energies E j are fixed in all analyses exceptthose of the background ROI data for constructing our lists of instrumental lines; in those analyses only, the E j are also modelparameters. would be to allow the GP model to have an energy-dependent mean. Equivalently, we could include a parametricmodel component (such as a power-law or exponential) to model the clear upward trend in the data at low energiesobserved in e.g. Fig. S1, with the GP model then describing fluctuations about that parametric component. Suchan approach would likely result in smaller values of the hyperparameter A GP and, potentially, increased sensitivity.Such a hybrid parametric plus GP modeling approach could be explored in future work.Our goal is to look for narrow lines on top of a smooth continuum flux. We know that even a narrow line willmanifest itself as a broader feature in the detector-level data due to the detector response. So the correlation-lengthof the GP kernel has a lower-bound set by the detector resolution. Because the energy resolution δE of XMM-Newton increases linearly with energy ( i.e. , δE/E is roughly constant), a stationary kernel with a fixed correlation lengthis not adequate and a kernel of the form in (S3) is more natural. However, we expect the continuum to be muchsmoother, even before the smearing induced by the detector resolution. A common approach in GP literature whenthe hyperparameters are not motivated from some other considerations is to fit them to the data. This approach leadsto the best-fit values σ E ≈ .
608 ( σ E ≈ .
77) for MOS (PN) in the first annulus, with comparable results in the annulifurther from the Galactic Center. However, we chose to fix σ E = 0 . σ E the GP model isable to capture smaller-scale fluctuations in the data, absorbing what would otherwise be attributed to narrow lines.Lastly, note that while the marginal likelihood in (S2) is defined within the context of Bayesian statistics, as it isobtained by integrating the likelihood times prior distribution for the formal GP model parameters, we will use thelikelihood to perform frequentist parameter inference. This approach is called the “hybrid approach” in [42]. As notedin [42], the asymptotic expectations for the distribution of the TS constructed from the marginal likelihood may differfrom the frequentist expectations [71], because of the use of the Bayesian marginal likelihood, and so in principlethe p -values and upper-limit criteria should be calibrated on Monte Carlo (MC). However, as we show below, wefind through MC simulations that in our examples the TS statistical distributions follow the asymptotic frequentistexpectations to high accuracy. With that in mind, we briefly review the asymptotic expectations for translatingdiscovery TS values to p -values and forming 95% one-sided upper limits.As discussed in the main Letter, the TS in favor of the signal model is given by t = − (cid:2) max θ log p ( d k | θ ) − max θ nuis log p ( d k |{ A sig = 0 , θ nuis } (cid:3) , (S4)where the second term is the maximum marginal likelihood for the null model without a signal line. When searchingfor evidence of DM, the discovery TS is set to zero for unphysical model parameters (in that case, sin (2 θ ) < χ distributed with one degree of freedom under the null hypothesis(see, e.g. , [71]). In addition to searching for evidence of the signal model over the null hypothesis using t , we also set95% one-sided upper limits using the likelihood ratio. We define the profile likelihood ratio q ( A sig ) by q ( A sig ) = − (cid:2) max θ nuis log p ( d k |{ A sig , θ nuis } ) − max θ log p ( d k | θ ) (cid:3) , (S5)where in the first term we maximize the log-likelihood over the nuisance parameters θ nuis at fixed signal parameter A sig . Let ˆ A sig be the best-fit signal parameter; i.e. , q ( ˆ A sig ) = 0. Then, the 95% one-sided upper limit A is givenby the value A > ˆ A sig which satisfies q ( A ) ≈ − .
71 [71].Note that the profile likelihood in (S5) is computed as a function of the signal normalization A sig at fixed UXLenergy (or, equivalently, fixed DM mass). All of the analyses presented in the main Letter are performed in thisway ( i.e. , we have a grid of different UXL energies to probe and then for each fixed energy we compute the profilelikelihood as a function of the signal-strength amplitude). In SM Sec. III A, however, we consider our ability to localizea putative signal in m a and sin (2 θ ) using synthetic data. In that analysis, and that analysis only, we simultaneouslyconstrain the mass and the signal strength.
2. Instrumental and astrophysical background lines
Several instrumental and astrophysical background spectral lines are expected to contribute to the observed fluxspectra. Here, we outline the procedure by which we obtain our candidate instrumental and astrophysical backgroundlines, with the complete list of included lines for MOS and PN presented in Tab. S1 and Tab. S2, respectively.We adopt an initial instrumental line list for PN and MOS from Refs. [72, 73]. We then analyze the stacked data,for MOS and PN independently, in the background ROI in order to test for the presence of each candidate line. Weuse an analysis framework analogous to that we use in the background-subtracted signal ROI data: in particular, our
Energy [keV] . . . . F l u x [ ph o t o n s / s / k e V ] MOS Annulus 1
Signal SpectrumBkg. SpectrumBkg-Sub. Spectrum
Energy [keV]
MOS Annulus 8
Energy [keV] . . . . . . F l u x [ ph o t o n s / s / k e V ] PN Annulus 1
Energy [keV]
PN Annulus 8
Figure S1. Examples of the signal region spectra for MOS (top panels) and PN (bottom panels) in Ring 1 (left panels) andRing 8 (right panels) with and without background subtraction in red and black, respectively. The background-region spectraare shown in grey. Many of the large instrumental features that are removed when looking at the background-subtracted data.Note that for visual clarity these spectra have been down-binned by a factor of 4. analysis of the background ROI data incorporates GP modeling for the continuum emission, in addition to the set ofputative instrumental lines. We test for known instrumental lines in the vicinity of: 4.51 keV (Ti Kα ), 5.41 keV (CrK α ), 5.90 keV (Mn K α ), 5.95 keV (Cr K β ), 6.40 keV (Fe K α ), 6.49 keV (Mn K β ), 7.06 keV (Fe K β ), 7.47 keV (NiK α ), 8.04 keV (Cu K α ). During this process, we allow the central location of the background lines to float by up to25 eV. Lines which are detected with t >
16 (4 σ local significance) in the background data analysis are accepted attheir best-fit energy as a new component of our residual background model. In MOS, we accept instrumental lines atenergies: 5.42 keV, 5.915 keV, 6.425 keV, 7.07 keV, 7.485 keV, and 8.06 keV. In PN, we accept instrumental lines at4.52 keV , 5.42 keV, 5.95 keV, and 6.39 keV.After constructing our list of instrumental background lines we include them in our analyses of the signal-ROIbackground-subtracted data sets. In particular, each line is given an intensity nuisance parameter in each ring.Given our procedure of subtracting the background flux from the signal region, variability in the instrumental linesbetween observations can result in the best fit instrumental line intensity in our data set having a positive or negativenormalization. Accordingly, we allow the normalization of the instrumental lines to be either positive or negative.We also develop an initial list of candidate astrophysical background lines following [8], by selecting those withemissivities greater than 5 × − photons/cm /s at a temperature of 1 keV, which is the approximate temperatureof the hot component of the Galactic Center, using the AtomDB database [74]. We include additional iron lines thatare known to produce emission in the inner Galaxy [75]. Taking this preliminary list, we then inspect the innermostring and determine all lines which appear with TS t > Energy[keV] Origin Type Ring 1 Ring 2 Ring 3 Ring 4 Ring 5 Ring 6 Ring 7 Ring 82.46 S Astro. χ associated with the addition/removal of the line fromthe best-fit background model which is obtained after our line-dropping procedure. Bolded values indicate the inclusion of aline in a ring’s background model.Line Energy[keV] Origin Type Ring 1 Ring 2 Ring 3 Ring 4 Ring 5 Ring 6 Ring 7 Ring 82.46 S Astro. χ associated with the addition/removal of the line fromthe best-fit background model which is obtained after our line-dropping procedure. Bolded values indicate the inclusion of aline in a ring’s background model. are detected with moderate significance (we use the criteria t >
3) in the background-subtracted data set. Note thatin Tabs. S1 and S2 we indicate whether the line is included in each annulus. To determine the significance of agiven line we proceed iteratively, starting with the full list of lines and then calculating the change in the maximumlikelihood when a given line is removed from the model.In Fig. S2 we illustrate example fits for our fiducial analyses to the data without the inclusion of an UXL. Thesefits are to the same background-subtracted data as illustrated in Fig. S1, as labeled. In black we show the combinedbest-fit model, which is the sum of the GP model contribution (dark red) and the contributions from the individualastrophysical and instrumental lines (colored curves). Note that the number of background lines differs between eachof the annuli because the important background lines are determined independently for each annulus.
3. The joint likelihood and background mismodeling
After constructing the profile log likelihoods q k ( A sig ) in each energy annulus ( k = 1 , , · · · ,
8) we then use (1)to convert from A sig , which has units of cts / cm / s / sr, to sin (2 θ ). To do so we use the background-subtracted D - − . . . . . . F l u x [ ph o t o n s / s / k e V ] MOS Annulus 1 − . . . . . . MOS Annulus 8
Energy [keV] . . . . F l u x [ ph o t o n s / s / k e V ] PN Annulus 1
Energy [keV] . . . . PN Annulus 8
Figure S2. The same background-subtracted data sets illustrated in Fig. S1 (also down-binned), but now shown along withtheir best-fits under the null hypothesis. The best-fit model prediction is shown in black, which may be decomposed into thecontribution from the GP model (dark red) and the contributions from the individual background lines (colored curves). Notethat the background lines to include in the analysis are determined independently in each annulus, as described in the text. factors, as discussed in the main Letter. Then, at each test mass point for the DM model we construct the jointprofile likelihood q joint (cid:0) sin (2 θ ) (cid:1) = (cid:88) k =1 q k (cid:0) sin (2 θ ) (cid:1) , (S6)for both MOS and PN independently. Later, we will also combine the MOS and PN profile likelihoods to constructour final joint profile likelihood that we use to search for evidence of decay DM. First, however, we analyze the jointMOS and PN profile likelihoods independently for evidence of background mismodeling.We test and account for possible background mismodeling by extending the background model to include a compo-nent that is totally degenerate with the signal. This is a conservative approach that would remove all sensitivity to aUXL if the amplitude for this additional signal-like component were left free. Therefore we penalize the amplitude ofsuch a signal like feature in the background model with a zero-mean Gaussian likelihood with variance hyperparameter σ . The approach we follow was developed and implemented in [76–78] within the context of searches for narrowspectral features in γ -ray astronomy and in the context of the Higgs boson search by the ATLAS experiment, whereit is called the “suprious signal” [40]. We extend the likelihood to include two “spurious signal” nuisance parameters,one for the MOS data and one for the PN data. The MOS and PN likelihoods are then combined to produce the jointlikelihood that we use for probing the DM model.After extending the background model to include a signal-like component constrained by σ , the resulting profilelikelihood (for either the MOS or PN data) is given by˜ q joint (cid:0) sin (2 θ ) (cid:1) = max A spur (cid:34) q joint (cid:0) sin (2 θ ) + A spur (cid:1) − ( A spur ) σ (cid:35) , (S7)where q joint is defined in (S6). Note that the profile likelihood now depends on the hyperparameter σ , whichdetermines the strength of the spurious-signal nuisance parameter For example, in the limit σ → A spur →
0) and the modified profile likelihood ˜ q approaches the un-modifiedlikelihood q . However, in the opposite limit σ → ∞ we completely lose constraining power and ˜ q joint (cid:0) sin (2 θ ) (cid:1) → (2 θ ).In practice, we determine the value of the hyperparameter at each test mass point independently. The philosophyis that if there is evidence that the background model is not properly describing the data in the immediate energyside-bands around a mass point of interest, then we should account for the possibility, through A spur , of similarbackground mismodeling at our mass point of interest. Specifically, we implement the following approach. At a givenmass point m mχ , where m is the index that labels the mass point, we consider the subset of test mass points in a 2 keVwindow around m mχ , masking: (i) a 0.4 keV window in mass around m mχ and (ii) masking 0.1 keV windows aroundthe locations are all background lines that were included in the analyses of the annuli. Each test mass point withinthis side-band window has a best-fit sin (2 θ ) from the likelihood analysis without the inclusion of the spurious-signalnuisance parameter. The ensemble of best-fit points in the side-band window is denoted by { sin (2 θ ) } m . We computethe variance over this ensemble of best-fit points, Var (cid:2) { sin (2 θ ) } m (cid:3) observed . The observed variance is then comparedto the expected variance Var (cid:2) { sin (2 θ ) } m (cid:3) expected , and specifically we set σ , m = max (cid:104) , Var (cid:2) { sin (2 θ ) } m (cid:3) observed − Var (cid:2) { sin (2 θ ) } m (cid:3) expected (cid:105) , (S8)where σ , m denotes the hyperparameter at the mass point m mχ . The expected side-band best-fit varianceVar (cid:2) { sin (2 θ ) } m (cid:3) expected is computed from 500 MC simulations of the null hypothesis. The null hypothesis model isthat given by the fit of the background model to the data without any extra UXL signal components.We expect σ , m to be non-zero if there is background mismodeling in the energy side-band, which increases thevariance of observed best-fit points relative to the expectation under the null hypothesis. However, sometimes σ , m will be non-zero simply because of statistical fluctuations in the observed side-bands, in which case the nuisanceparameter will weaken the limits more than intended. However, this occasional weakening of the limits is worth theadvantage of having an analysis framework that is more robust to mismodeling. Indeed, we know that there is anopportunity for some degree of background mismodeling because we have chosen to only include background lines thatpass some significance threshold, and thus the aggregate effect of the sub-threshold lines could lead to mismodelingthat could be partially mitigated by A spur . Energy [keV] − − − − σ s i n ( θ ) MOS Stat.MOS Sys. PN Stat.PN Sys
Figure S3. The spurious-signal hyperparameter σ , m (labeled MOS Sys. and PN Sys.), as computed in (S8), as a functionof the UXL energy. For both MOS and PN the nuisance parameter A spur is predominantly active at low energies, and it playsa more significant role in PN than in MOS. We compare the hyperparameter to the statistical uncertainties (labeled MOSStat. and PN Stat.), which are computed from the Hessian of the log-likelihood (without the spurious-signal) about the best-fitmixing angle at a fixed energy. In Fig. S3 we illustrate the values of σ , m (labeled MOS Syst. and PN Syst.) that we find from the data analysesof the MOS and PN data. We compare the hyperparameter to the statistical uncertainty on sin (2 θ ), labeled MOSStat. and PN Stat. Note that the statistical uncertainties are computed from the Hessian of the log-likelihood, forthat data set, about the best-fit coupling at a fixed UXL energy, without the inclusion of A spur . For both MOSand PN we see that the background mismodeling uncertainties, as captured by A spur , may dominate the statisticaluncertainties at some low energies, though the nuisance parameter appears more important for PN than for MOS.It is interesting to consider the ensemble of discovery TSs in favor of the DM model across all tested mass points.We denote this distribution of TSs without the spurious-signal by T , while with the inclusion of the spurious-signalnuisance parameter we call this distribution T sys . We expect T sys to have fewer high-TS points than T . In the leftand center panels of Fig. S4 we illustrate the distributions of TSs for both T (labeled Data) and T sys (labeled Dataw/ Nuisance Parameter) for MOS and PN, respectively. More specifically, in that figure we illustrate the survivalfractions for the distributions, which show the faction of TSs in T or T sys with a value above the TS indicated on the x -axis. Asymptotically we expect, up to the caveat that we used the Bayesian marginal likelihood of the GP to defineour TSs, that the TSs should be χ distributed [71]. The survival function of the χ distribution is shown in Fig. S4.We verify with a large number of MC simulations that the that the null-distribution of TSs is indeed chi -distributedfor both MOS and PN datasets. The results of these tests are labeled “Monte Carlo” in Fig. S4 and overlap with the χ distribution, providing evidence that we are in the asymptotic regime [71]. − − t − − − Su r v i v a l F r a c t i o n MOS
DataData w/ NuisanceMonte Carlo χ Distribution − − t PN 10 − − t Joint
Figure S4. (Left) The survival function of the test statistic for discovery in the analysis of the MOS data. Under the nullhypothesis, and for a large number of samples, the survival fractions are expected to follow the χ distribution, as verified byMC (as labeled). At a finite number of samples the expected chi-square distributions are found from MC to be expected tobe contained within the green and gold shaded regions at 68% and 95% confidence, respectively. The negligible effect of thesystematic nuisance parameter can be seen by comparing the survival function without the nuisance parameter (red, labelled“Data”) and with the nuisance parameter (blue, labeled “Data w/ Nuisance Parameter”). (Center) As in the left panel, but forthe PN analysis. (Right) The survival function for the joint analysis of MOS and PN data. In blue, the survival function forthe joined PN and MOS analysis without systematic nuisance parameters; in red, the survival function for the joint analysiswhen the PN and MOS results are corrected by their independently-tuned systematic nuisance parameters prior to joining. Because there are a finite number of samples in T and locations spaced within the detector energy resolution arecorrelated, the survival function for the observed data is not expected to follow the χ -distribution exactly. Thegreen and gold bands in Fig. S4 show the 68% and 95% containment regions for the survival fraction computed over500 MC realizations of T . We expect that the data should fall within these bands if no signal is present, which isanalogous to the green and gold bands for the significance in Fig. S5. In the left and center panels of Fig. S4 we maysee that the distributions of T for MOS and PN are broadly consistent with the MC expectations. The distributionsof T sys , as expected, fall off slightly faster at large values of the TS. The right panel of Fig. S4 shows the survivalfraction for the combined analysis where we combine the MOS and PN profile likelihoods, with and without thespurious-signal. The most significant test point has a significance slightly above 2 σ local significance, which is lessthan 1 σ in global significance. Thus, we conclude that there is no evidence for decaying DM in our analysis aboveour 5 σ global predetermined detection threshold.The effect of the spurious-signal nuisance parameter on the individual MOS and PN limits is illustrated in Fig. S5.The inclusion of the nuisance parameter slightly decreases the discovery TSs at low masses and also causes a slightweakening of the limits. Note that the expectations under the null hypothesis are indicated for the spurious-signal-corrected analysis in that figure. II. EXTENDED RESULTS
In this section we present extended results for the analyses that go into producing Fig. 4. First, we provide thebest-fit normalizations for our GP kernels, presented in Tab. S3. In Fig. S6 we present the main result in Fig. 4 interms of the DM lifetime instead of in terms of sin (2 θ ). The result in Fig. S6 is more general than in Fig. 4 since itholds for more general DM models beyond the sterile neutrino model. Note, however, that this figure applies to DM − − − − s i n ( θ ) DM under production (4)(3)(2)(1) DM under production (4)(3)(2)(1) m χ [keV] − S i g n i fi c a n ce ( √ t ) w/ syst. nuis. param. w/o syst. nuis. param. m χ [keV] Figure S5. As in Fig. 4, but for the MOS (left panel) and PN (right panel) analyses individually and with and without thespurious-signal nuisance parameter. Note that the 1 σ and 2 σ expectations for the limits are shown for the case with thespurious-signal nuisance parameter. The limits without the nuisance parameter are slightly stronger at low masses. m χ [keV] τ [ s ] (4)(3)(2)(1) Figure S6. As in Fig. 4, but interpreted as limits on the DM lifetime. This figure applies for DM whose decays produce a singlemono-energetic photon at energy m χ /
2. If the DM decay produces two photons (as in an axion model), then the lifetime limitsare twice as strong. whose decays produces one mono-energetic photon at energy m χ /
2. Axion-like models produce two photons duringthe decay, in which case the limits are twice as strong as those shown in Fig. S6.
Instrument Ring 1 Ring 2 Ring 3 Ring 4 Ring 5 Ring 6 Ring 7 Ring 8MOS 6 . × − . × − . × − . × − . × − . × − . × − . × − PN 2 . × − . × − . × − . × − . × − . × − . × − . × − Table S3. The best fit normalization of the GP kernel for each ring in PN and MOS. We present √ A GP in units ofphotons/cm /s/keV for A GP defined in (S3). Next, we show our results from the analyses to the individual MOS and PN annuli. In Fig. 3 we illustrate thebest-fit flux and significance (times the sign of the excess or deficit) for the UXLs for the innermost MOS annulus. InFigs. S7 through S14 we show the corresponding figures for the rest of the annuli and for PN. Note that the shadedgrey regions denote the masks that we use to avoid searching for UXLs in the direct vicinity of background linesincluded in the analyses.The distribution of discovery TSs that we observe in the individual annuli all appear consistent with expectationsfrom MC, as illustrated in Fig. S15 for MOS and Fig. S16 for PN. These figures illustrate the survival fractions of TSs,0 − . . . F l u x [ c t s / c m / s / s r ] MOS Annulus 1 E [keV] − S i g n i fi c a n ce ( √ t ) PN Annulus 1 E [keV] Figure S7. The left panel is a reproduction of Fig. 3, while the right panel illustrates the same quantities for the fit to theinnermost PN annulus. − . . . F l u x [ c t s / c m / s / s r ] MOS Annulus 2 E [keV] − S i g n i fi c a n ce ( √ t ) PN Annulus 2 E [keV] Figure S8. As in Fig. S7 but for annulus 2. as in Fig. S4, but at the level of the individual annuli instead of the joint analysis. Note that the MC expectations areconstructed independently for each annulus and each data set. These results do not include the systematic nuisanceparameter since that is only included at the level of the joint likelihood, after combining the results from all of theindividual annuli.1 − . . . F l u x [ c t s / c m / s / s r ] MOS Annulus 3 E [keV] − S i g n i fi c a n ce ( √ t ) PN Annulus 3 E [keV] Figure S9. As in Fig. S7 but for annulus 3. − . . . F l u x [ c t s / c m / s / s r ] MOS Annulus 4 E [keV] − S i g n i fi c a n ce ( √ t ) PN Annulus 4 E [keV] Figure S10. As in Fig. S7 but for annulus 4.
III. SYNTHETIC SIGNAL TESTS FOR THE FIDUCIAL ANALYSIS
In this section we verify that our analysis framework has the ability to discover real signals if they are present inthe data. We do so by injecting a synthetic signal into the real data and analyzing the hybrid data set with ourfull analysis. We also demonstrate the full analysis as applied to fully synthetic data generated with varying injectedsignal strengths.2 − . . . F l u x [ c t s / c m / s / s r ] MOS Annulus 5 E [keV] − S i g n i fi c a n ce ( √ t ) PN Annulus 5 E [keV] Figure S11. As in Fig. S7 but for annulus 5. − . . . F l u x [ c t s / c m / s / s r ] MOS Annulus 6 E [keV] − S i g n i fi c a n ce ( √ t ) PN Annulus 6 E [keV] Figure S12. As in Fig. S7 but for annulus 6.
A. Signal injection in real data
For injection tests in the real data, we chose a DM mass m χ = 7 . (2 θ ) = 2 . × − .We chose this mixing angle because we expect such a signal to be detected at approximately 5 σ significance. Weforward model this signal through the appropriate MOS and PN detector responses, draw Poisson counts, and thenadd these counts to the actual data sets. The results of the data analysis of the hybrid data are illustrated in Fig. S17.In the top panel we show the 95% upper limits for MOS, PN, and the joint analysis, with and without the systematicnuisance parameter. Note that the injected signal is indicated by the red star. The upper limits weaken at theinjected signal point, as expected, and do not exclude the injected signal coupling. In the second row we show thecorresponding detection significances. The signal is detected at nearly 5 σ in MOS alone and at around 2 σ in PN.3 − . . . F l u x [ c t s / c m / s / s r ] MOS Annulus 7 E [keV] − S i g n i fi c a n ce ( √ t ) PN Annulus 7 E [keV] Figure S13. As in Fig. S7 but for annulus 7. − . . . F l u x [ c t s / c m / s / s r ] MOS Annulus 8 E [keV] − S i g n i fi c a n ce ( √ t ) PN Annulus 8 E [keV] Figure S14. As in Fig. S7 but for annulus 8.
The systematics nuisance parameter slightly reduces the significance of the discovery, but by a minimal amount sincewe mask a 0.4 keV window around the test mass when determining the systematics nuisance parameter. In the thirdrow we show how the discovery of the injected signal extends the survival function to higher TS values. Lastly, inthe bottom row we show the 1, 2, and 3 σ best-fit regions in the m χ -sin (2 θ ) plane for the DM candidate. In red wemark the location of the injected signal, which is recovered appropriately.4 − − − Su r v i v a l F r a c t i o n MOS Annulus: 1DataMonte Carlo χ Distribution MOS Annulus: 2 MOS Annulus: 3 MOS Annulus: 4 − − t − − − Su r v i v a l F r a c t i o n MOS Annulus: 5 − − t MOS Annulus: 6 − − t MOS Annulus: 7 − − t MOS Annulus: 8
Figure S15. As in Fig. S4 but for the individual MOS annuli. Note that the systematic nuisance parameter has not beenapplied since that is only incorporated in the joint likelihood that combines the results from the individual annuli. − − Su r v i v a l F r a c t i o n PN Annulus: 1DataMonte Carlo χ Distribution PN Annulus: 2 PN Annulus: 3 PN Annulus: 4 − − t − − Su r v i v a l F r a c t i o n PN Annulus: 5 − − t PN Annulus: 6 − − t PN Annulus: 7 − − t PN Annulus: 8
Figure S16. As in Fig. S15 but for the PN data sets.
B. Signal injection in synthetic background data
For injection tests on the real data, we first generate synthetic data according to the best-fit null models for eachof the eight rings studied in MOS and PN data sets. We then inject a synthetic signal at a specified value of sin (2 θ )on top of the null-model realizations using the same procedure as applied for the signal injection on the real data andrepeat our full analysis procedure in search of the injected signal with the exception that we do not apply a nuisanceparameter tuning and correction. We perform 1000 independent realizations and analyses for each value of sin (2 θ ,and we repeat this procedure for 30 values of sin (2 θ ) between 10 − and 10 − for two different neutrino masses: 7.05 m χ [keV] − − − − s i n ( θ ) (4)(3)(2)(1) MOS
Limit w/ NuisanceLimit w/o Nuisance m χ [keV] (4)(3)(2)(1) PN m χ [keV] (4)(3)(2)(1) Joint m χ [keV] − . − . . . . S i g n i fi c a n ce ( √ t ) w/ Nuisancew/o Nuisance m χ [keV] m χ [keV] − − t − − − Su r v i v a l F r a c t i o n DataData w/ NuisanceMonte Carlo χ Distribution − − t − − t . . . . m χ [keV] − × s i n ( θ ) Injected SignalRecovered Signal . . . . m χ [keV] . . . . m χ [keV] Figure S17. The results of the analysis of the hybrid data that consists of the real MOS and PN data plus a synthetic DMsignal. The DM signal is generated with mass m χ = 7 . (2 θ ) = 2 . × − as described in the text.The top, middle, and third rows are analogous to Figs. S4 and S5, but for the hybrid data set. The last row shows the 1, 2, and3 σ recovered parameter space for the signal in the mass and mixing angle plane. The best-fit recovered signal is indicated indark blue, while the red star denotes the true value injected. The synthetic signal is appropriately recovered, adding confidencethat our analysis procedure has the ability to detect real DM signals if present in the data. keV and 11.5 keV. The results of the data analysis of the hybrid data are illustrated in Fig. S18. In the top row, weshow ensemble statistics for the 95% upper limits as a function of injected signal strength for the two neutrino massesstudied in this test. In the bottom row, we show the ensemble statistics of the recovered detection test statistic asa function of the injected signal strength. The upper limits weaken with increasing injected signal strength withoutexcluding the true value of the injected signal. Moreover, the detection test statistic smoothly increases as a functionof increasing injected signal strength. Critically, at large injected signal strength, the test statistic safely exceeds T S ≈
30, which is the approximate threshold for a 5 σ detection after correcting for the look-elsewhere effect.6 − − s i n ( θ r ec ) m χ = 7 keV Median Limit1 / σ Containment m χ = 11 . − − sin (2 θ inj )10 − − t Median Discovery TS1 / σ Containment − − sin (2 θ inj ) Figure S18. (
Top Row ) In red, the median 95 th percentile upper limit on the recovered signal as a function of the injection signalstrength at two neutrino masses evaluated on synthetic data. We additionally indicate the 1 and 2 σ containment intervals forthe ensemble of upper limits realized at each injected signal strength. Note that these upper limits are not power constrained.These results demonstrate that our analysis framework places robust upper limits that do not rule out an injected signal.( Bottom Row ) In black, the median recovered detection test statistic for a signal injected in the synthetic data as a function ofthe injected signal strength, with the 1 and 2 σ containment intervals also indicated. Under the null hypothesis, the detectiontest statistic should follow a χ -distribution; the median and 1 σ and 2 σ percentile values of the χ -distribution are indicatedby dashed grey lines. These results demonstrate that our detection test statistic follows its theoretically expected distributionunder the null hypothesis (sin (2 θ inj ) = 0) and that our analysis framework can robustly identify a signal which is present inthe data. The results are smoothed with a Savitzky–Golay filter for clarity. IV. SYSTEMATIC ANALYSIS VARIATIONS
Our fiducial result, which is illustrated in Fig. 4, made a number of physics-level and analysis-level choices. Thesechoices are justified in the main text and the supplementary results of the proceeding sections of the SM. Still, it isworthwhile to consider how our results change for different physics and analysis assumptions and choices, as this givesan indication of the robustness of the limits and significances shown in Fig. 4.
A. Alternate DM Density Profiles
In the main text, and in particular in Fig. 4, we adopted the conservative DM profile that was shown in Fig. 1.As already described, the present expectation is that in the absence of baryons, the DM halo is well described by anNFW profile. Baryons are then expected to contract this profile, increasing the DM density towards the GC, andpotentially also introducing a core on top of this. For our fiducial analysis we conservatively assumed an uncontractedNFW halo, using the most conservative parameters determined within the 68% best fit region of Cautun et al. [57].In particular, we used an NFW profile with r s = 19 . ρ DM = 0 . − − − − s i n ( θ ) DM under production (4)(3)(2)(1) m χ [keV] − S i g n i fi c a n ce ( √ t ) Fiducial. Cautun best fit Contracted
Figure S19. As in Fig. 4, but for three different DM density profiles, all based upon Ref. [57]. In solid curve we show ourfiducial results, corresponding to the uncontracted NFW profile with a conservative density. The dashed curve then shows ourresults using the best fit NFW profile, whereas in dashed we show the stronger limits that would be obtained with a contractedDM distribution. Details of the distributions are provided in the text.
GeV/cm .In Fig. S19, we show our main results if instead we repeat the analysis for the best fit NFW profile determinedin Ref. [57], which corresponds to r s = 15 . ρ DM = 0 .
31 GeV/cm , as well as showing results for the morerealistic contracted profile. There is not a parametric form for the contracted profile, however, Ref. [57] providesa best fit model for the DM mass distribution, which we use to infer the density and then D -factor. The modelonly provides an estimate down to 1 kpc from the GC, within which we conservatively assume the density profile iscompletely cored.As the figure demonstrates, adopting a more realistic contracted DM profile strengthens our limits by roughly afactor of 2. Importantly, however, changing the profile does not appreciably change the distribution of significance,and we continue to see no clear evidence for an UXL. B. Dependence on the GP model
For our fiducial analysis we use the GP kernel given in (S3) with the choice σ E = 0 .
3. This choice was made sothat the residual background model has the ability to adjust on scales around one order of magnitude larger in scalethan the energy resolution of the detectors, which are δE/E ∼ .
03. In this section we verify that our results do notdepend in detail upon the particular value chosen.First, we consider a small modification to our default analysis by taking σ E = 0 . σ E = 0 .
4. The resultsof these analyses are shown in Figs. S20 and S21. In those figures we show the 95% upper limits (upper panel),significances (middle panel), and survival fraction of significances (bottom panel). We give the results both with anwithout nuisance parameters. There is a slight trend where increasing σ E leads to a corresponding strengtheningof the sensitivity, though this difference is minor compared to other choices in the analysis. In general, the resultsappear robust to the choice of σ E .Next, we consider changing the GP modeling more significantly by adopting an alternate kernel. In particular, weconsider the standard (and stationary) double exponential kernel K ( E, E (cid:48) ) = A GP exp (cid:20) − ( E − E (cid:48) ) σ (cid:21) , (S9)which has the hyperparameter σ . Note that our fiducial kernel, given in (S3), has the property whereby the correlation8 m χ [keV] − − − − s i n ( θ ) (4)(3)(2)(1) MOS
Limit w/ NuisanceLimit w/o Nuisance m χ [keV] (4)(3)(2)(1) PN m χ [keV] (4)(3)(2)(1) Joint m χ [keV] − S i g n i fi c a n ce ( √ t ) w/ Nuisancew/o Nuisance m χ [keV] m χ [keV] − − t − − − Su r v i v a l F r a c t i o n DataData w/ NuisanceMonte Carlo χ Distribution − − t − − t Figure S20. The analogues of Figs. S4 and S5, but changing the kernel correlation length to σ E = 0 . σ E = 0 . σ E = 0 . . length increases with the energy resolution of the detector. The kernel in (S9), on the other hand, has a fixed correlationlength as a function of energy. In Figs. S22 and S23 we show the results of using the double exponential kernel withscale length σ = 0 . σ = 1 . σ slightly increases the limits. However, the differences between the double-exponential kernel results andour fiducial results are minor and most evident at high DM masses, m χ , where the two kernels predict the largestdifferences. In particular, we find no evidence for decaying DM with the alternate kernels and similar 95% upperlimits. The systematic uncertainty associated with this choice is generally less than other aspects of the analysis suchas our assumptions regarding the DM density profile. C. Analysis of Fully Stacked Data
In the main body, we divided our signal ROI into rings and modeled the flux independently in each ring. Themotivation behind this choice was to incorporate spatial information into the analysis, particularly as we expect theflux of an actual DM decay signal to steadily increase towards the GC. Here we show the results of an alternativeapproach where instead of modeling the data ring-by-ring, we instead combine the data in the innermost three ringsof the signal ROI and model that directly. We effectively are then left with a single combined ring, which we analyzeusing our fiducial procedure.In Fig. S24 we show the resulting limit in the case where we also subtract the background-ROI flux from the stackedsignal region data. While there are small differences, the resulting sensitivity and limits from this simpler approachare in good qualitative agreement to those of our default analysis. In detail, the result here are slightly weaker, whichis as expected because there is less information in the signal ROI (we use fewer rings and by stacking the spatial9 m χ [keV] − − − − s i n ( θ ) (4)(3)(2)(1) MOS
Limit w/ NuisanceLimit w/o Nuisance m χ [keV] (4)(3)(2)(1) PN m χ [keV] (4)(3)(2)(1) Joint m χ [keV] − S i g n i fi c a n ce ( √ t ) w/ Nuisancew/o Nuisance m χ [keV] m χ [keV] − − t − − − Su r v i v a l F r a c t i o n DataData w/ NuisanceMonte Carlo χ Distribution − − t − − t Figure S21. As in Fig. S20 but with σ E = 0 .
4. The limit is slightly strengthened, although again the differences are notsignificant. information is partially erased).Next, in Fig. S25 we repeat this procedure but without subtracting the background flux. The differences arenow more noticeable - the expected and resulting limit undergoes larger fluctuations and there are several mildlysignificant excesses. This emphasizes the importance of the background subtraction procedure in simplifying thedata, particularly around bright instrumental lines.
D. Parametric Modeling without Background Subtraction
In this section, we detail an alternate analysis to the one presented in the main body of this paper and providea comparison between the fiducial and alternate analysis in a representative example over the 8-9 keV mass range.The alternate framework uses the same data as used in our primary analysis. However, a more traditional approachis adopted for the background modeling. Firstly, we consider the unsubstracted data in each ring within the signalROI. The flux within each annulus is then modeled as follows. The background and putative signal lines are treatedidentically to our fiducial approach, but the smooth background contribution is modeled parametrically using anunfolded second order polynomial, rather than with a GP model. The three parameters that define the quadraticbackground component are treated as nuisance parameters and profiled over. As the quadratic background has lessfreedom than the GP model, we restrict to a smaller energy range. Specifically, we determine the energy range byfitting a Gaussian to the detector response at a given putative signal energy, and we define our energy range to extend5 standard deviations out from the signal energy in either direction. In the 8-9 keV DM mass range, this correspondsto an approximate energy range of 0.60 keV. Furthermore, background lines within 7 standard deviations of the signalenergy are included in the model. Thus, in the 8-9 keV mass range, the only line included is the 4.52 keV instrumental0 m χ [keV] − − − − s i n ( θ ) (4)(3)(2)(1) MOS
Limit w/ NuisanceLimit w/o Nuisance m χ [keV] (4)(3)(2)(1) PN m χ [keV] (4)(3)(2)(1) Joint m χ [keV] − S i g n i fi c a n ce ( √ t ) w/ Nuisancew/o Nuisance m χ [keV] m χ [keV] − − t − − − Su r v i v a l F r a c t i o n DataData w/ NuisanceMonte Carlo χ Distribution − − t − − t Figure S22. As in Fig. S20 but with the alternate GP kernel, in (S9), with σ = 0 . line for all PN annuli. We do not include the systematic nuisance parameter modeling for this example.While our background modeling is significantly different in this case, we find again that our results are qualitativelyunchanged compared to our fiducial analysis. To provide a representative example, in Fig. S26 we show the comparisonbetween our fiducial analysis (without the systematics nuisance parameter to facilitate the comparison) and thisalternate approach over the mass range 8-9 keV. As can be seen, the expected sensitivity of the two approachesis almost identical. This is a significant further demonstration that our specific choice of background model is notunderpinning our sensitivity.1 m χ [keV] − − − − s i n ( θ ) (4)(3)(2)(1) MOS
Limit w/ NuisanceLimit w/o Nuisance m χ [keV] (4)(3)(2)(1) PN m χ [keV] (4)(3)(2)(1) Joint m χ [keV] − S i g n i fi c a n ce ( √ t ) w/ Nuisancew/o Nuisance m χ [keV] m χ [keV] − − t − − − Su r v i v a l F r a c t i o n DataData w/ NuisanceMonte Carlo χ Distribution − − t − − t Figure S23. As in Fig. S22 but with σ = 1 . m χ [keV] − − − − s i n ( θ ) (4)(3)(2)(1) MOS
Limit w/ NuisanceLimit w/o Nuisance m χ [keV] (4)(3)(2)(1) PN m χ [keV] (4)(3)(2)(1) Joint m χ [keV] − S i g n i fi c a n ce ( √ t ) w/ Nuisancew/o Nuisance m χ [keV] m χ [keV] − − t − − − Su r v i v a l F r a c t i o n DataData w/ NuisanceMonte Carlo χ Distribution − − t − − t Figure S24. The same results as presented in Figs. S4 and S5, however on a modified data set where instead of analyzing thesignal ROI divided into eight individual rings, we stack the inner three rings into a single annulus. As in our primary approach,we subtract the background ROI flux from the signal-region data. The results are comparable to, although slightly weakerthan, those from our fiducial approach, consistent with the reduced information available. m χ [keV] − − − − s i n ( θ ) (4)(3)(2)(1) MOS
Limit w/ NuisanceLimit w/o Nuisance m χ [keV] (4)(3)(2)(1) PN m χ [keV] (4)(3)(2)(1) Joint m χ [keV] − S i g n i fi c a n ce ( √ t ) w/ Nuisancew/o Nuisance m χ [keV] m χ [keV] − − t − − − Su r v i v a l F r a c t i o n DataData w/ NuisanceMonte Carlo χ Distribution − − t − − t Figure S25. As in Fig. S24, however considering the stacked signal ROI without subtracting the background. The limit isnoticeably worse, and several excesses appear, highlighting the importance of the background subtraction. − − − − s i n ( θ ) FiducialFiducial Expected AlternateAlternate Expected . . . . . . m χ [keV] − S i g n i fi c a n ce ( √ t ))