Surface brightness profile of the 3.5 keV line in the Milky Way halo
A. Boyarsky D. Iakubovskyi O. Ruchayskiy D. Savchenko




Lorentz Institute, Leiden University, Niels Bohrweg 2, Leiden, NL-2333 CA, The Netherlands

Discovery Center, Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, Copenhagen, DK-2100, Denmark
Bogolyubov Institute of Theoretical Physics, Metrologichna Str. 14-b, 03143, Kyiv, Ukraine

Discovery Center, Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, Copenhagen, DK-2100, Denmark

Bogolyubov Institute of Theoretical Physics, Metrologichna Str. 14-b, 03143, Kyiv, Ukraine



    We report a detection of 3.5 keV line in the Milky Way in 5 regions offset from the Galactic Center by distances from  10' to 35^∘. We build an angular profile of this line and compare it with profiles of several astrophysical lines detected in the same observations. We compare our results with other detections and bounds previously obtained using observations of the Milky Way. 


Surface brightness profile of the 3.5 keV line in the Milky Way halo
    D. Savchenko
    Received: date / Accepted: date
====================================================================



Dark matter (DM) is a universal substance, dominating mass of all galaxies and galaxy clusters and driving cosmological expansion of the Universe throughout most of its history.
It has evaded so far a definitive detection, although a number of astrophysical signals have been interpreted as originating from annihilations or decays of dark matter particles. 

Recently, an unidentified feature at E ∼ 3.5 keV in the X-ray spectra of dark matter dominated objects has been reported <cit.> (for review of the status of the line and its possible interpretations, see <cit.> and refs. therein).
Such a line could be interpreted as  coming from the decay of a DM particle with the mass ∼ 7 keV.
In most of these observations the detected line is very weak (its flux is usually 1-2% of the background) and its instrumental width exceeds the intrinsic one by several orders of magnitude. 
Therefore, it is not possible to check the origin of the line using the data from one object only, as one can always find an alternative explanation (statistical fluctuation, unknown astrophysical line; unmodelled instrumental feature).
The only possible approach is to study the morphology of the line, to compare its flux between different objects, and to study its surface brightness profile in each object. 
Of course, this approach is also subject to uncertainties and has to be applied with care as 
the DM content of astrophysical objects has not been determined with sufficient precision and DM distribution within each object even more so.
Nevertheless gathering more morphological data about the candidate line can help to do non-trivial and trustworthy consistency checks.



Line from the Milky Way.
A very important target for such morphological studies is our Galaxy.
The 3.5 keV line has been detected with
 observations inside the 30 pc (14') radius circle around the Galactic Center <cit.> (see also <cit.>). 
In all these papers, the line-like excess at 3.5 keV has been reported at high statistical significance. 
The strength of this line is consistent both with detection from the center of Andromeda galaxy and with non-detection from the blank sky dataset, where an upper limit on the flux at the level [0.7 × 10^-6]cts/sec/cm^2 has been reported <cit.>.
Alternative interpretations included K xviii  emission line complex at 3.47-3.52 keV <cit.> or Charge Exchange between neutral Hydrogen and bare sulphur ions leading to the enhanced transitions between highly excited states and ground states of S xvi <cit.>.
The Potassium interpretation is probably at tension with the recent observation of Perseus galaxy cluster by the Hitomi spectrometer <cit.> that did not reveal Potassium complex of lines at E∼ 3.5 keV.
These data is consistent with a broader DM line or a complex of S xvi transitions.

Ref. <cit.> has reported a detection of a 3.5 keV line from the Chandra Deep Field South and COSMOS fields (∼ 120^∘ off-GC) with NuSTAR telescope, while <cit.>  detected the 3.5 keV line from the same sky regions with the consistent flux via observations. 

The nature of 3.5 keV excess will probably be revealed with the launch of XRISM telescope <cit.>, that will be a true game changer, due to its superb spectral resolution. 
In the meanwhile we continue to investigate the surface brightness profile of the signal in the Milky Way using existing archival observations.

Our dataset.
We use 795 pointings of the  in the direction up to 35^∘ away from the GC (see the list of observations at <cit.>).
The total exposure exceeds 31.4 Msec (MOS cameras) and 8.9 Msec (PN camera).
We split observations into the regions (Reg1–Reg5) (see Table <ref>) and stack MOS and PN spectra for each of the spatial regions (using the same strategy we used in the previous works).
The resulting spectra are fit with a combination of astrophysical powerlaw continuum; unfolded powerlaw, representing cosmic-ray-induced background; and a set of narrow Gaussian lines, representing most prominent astrophysical and instrumental lines (the details of the data analysis are provided in the Supplementary material, Section <ref>).
The combined spectra for different regions, together with the best-fit models, are shown in Figs <ref>–<ref>.





Results.
 In all 5 spatial regions we have detected a 3.5 keV line. 
 Compared to our work <cit.> where the signal was extracted from the central 14', we have excluded the central 10' around the GC.
Nevertheless the signal was found in the annulus 10'-14' (Reg1) with the flux consistent with <cit.>.




Our results are summarized in Table <ref> and in Fig. <ref>. 
The local statistical significance of the line detection in the regions Reg1–Reg5 is ranging from 5.0σ (Reg5) to 2.1σ (Reg2).
The positions of lines in all 4 regions are consistent with the weighted average being
at E=[(3.49± 0.05)]keV.
The lines are found in both MOS and PN camera and their parameters are consistent. 
Therefore, in all regions we perform the combined fit to both cameras with the positions and flux per solid angle being the same. 
Along with the 3.5 keV line we have detected several other lines at the interval 3-4 keV: Ar XVII complex at ∼ 3.1 keV, Ca XIX complex at ∼ 3.9 keV and two weaker lines at 3.7 keV (possibly Ar XVII complex and/or instrumental Ca Kα line) and a line at 3.3 keV that is a combination of  Ar XVIII, S XVI and of K Kα instrumental line (c.f. <cit.>).

The angular profile of the  3.5 keV  line is shown in Fig. <ref>.
Having in mind possible DM interpretation of the signal, we would like to compare its distribution with a DM density profile of the Milky Way. 
The latter is not known with any reasonable precision in the inner region of our Galaxy. For example, Navarro-Frenk-White  models of the Milky Way halo typically have r_s > 100^∘, therefore our regions 1-5 are deeply inside r_s. It is well known that at such small distances from the center simple density profile models may not describe well even DM-only simulated halos, saying nothing about the real Galaxy.







Having this uncertainty in mind, we performed a combined fit to all 5 spatial regions with  relative normalization of the line in different region fixed in accordance with a Milky Way DM density profile. 
We explored several DM distributions in the Milky way (as described in Table <ref>). 
In each case we are able to find a good fit to the data with a significant (about 7σ) improvement for the quality of fit when adding a line at position E∼ 3.494 keV (see Table <ref> for details). 



The flux predicted by these profiles is also higher than the upper limit from the blank sky <cit.>.
Finally, we  compared the radial behavior of the 3.5 keV line  with that of strong astrophysical lines (see Fig. <ref>). We see that the flux in 3.5 keV line drops with the distance from the Galactic center slower than the flux in the astrophysical lines.



Discussion.
In this work we have detected the 3.5 keV line 
in 5 spatial regions with offsets from the GC ranging from 10' to 35^∘. We  tried to  fix the line flux ratios between different regions according to several DM density profiles (NFW, Einasto, Burkert) and were able to detect the line for every profile. The  best fitting profile among all that we tried is given by the NFW with r_s = [20]kpc <cit.> (see Table <ref>).
We also detected the line in all five regions independently (see Fig. <ref> and Figs. <ref>–<ref>). 
The fluxes in four regions (Reg2–Reg5)  are consistent with this best-fit NFW profile.
The upper bound reported in <cit.> is also consistent with the same DM density distribution. The best fit emission in the central part is somewhat higher than the prediction based on the other regions
(consistent with the results of <cit.>).
It is clear that DM density distribution at so small distances from the center (10 times smaller than r_s) is likely to be very different from NFW or any other simple model. It is is also possible that other potential sources of the 3.5 keV line (such as Potassium or Sulfur) contribute more significantly in the central region.

The prediction based on our best fit NFW model (as well as the upper bound from B14) is lower than the flux in the 3.5 keV line detected by the NuSTAR <cit.> and Chandra <cit.> satellites  in the COSMOS and Chandra Deep Field South (CDFS) fields. 
This disagreement requires further analysis. It could be related to cross-instrument
calibration or to underestimated errors on the flux. 

It is clear that the normalization of DM density profile (fixed e.g. by the local DM density around the Earth) is degenerate with the intrinsic decay width of DM particles (e.g. expressed in terms of the active-sterile neutrino mixing angle,  sin^2(2θ) for the case of sterile neutrino DM).
For the local DM density in the range [0.3-0.4]GeV/cm^3 <cit.> the sterile neutrino mixing angle sin^2(2θ) lies in the range (1.6 - 2.1) × 10^-11 fully consistent with the lower end of previous results <cit.>.

Our results are clearly not consistent with the bounds claimed in Ref. <cit.> using blank sky observations from the regions that partially overlap with our analysis. 
Although flux limits are not provided in <cit.>, we can convert the limit on the sterile neutrino mixing angle sin^2(2θ) into the upper bound on the flux.
Using NFW DM density profile used in <cit.>,
we estimate the upper limit on the flux to be

    F ≲ 3× 10^-3 cts/sec/cm^2/sr at 95% CL.

This would exclude the lines that we quote in Table <ref>.

In our view, this discrepancy is related to a very different approach to the signal and background modeling adopted in <cit.>. The authors of <cit.> limit the energy interval and  search for the line against a simple continuum model using the profile likelihood procedure.[Such a procedure is used, e.g., in searches for γ-ray lines by the Fermi-LAT collaboration <cit.>.]

The method used by <cit.> can be roughly compared to the sliding window procedure, where one subtracts from photon counts in a central bin Δ E counts from two adjacent bins with the width Δ E/2 
The sliding window method is based on the assumption that the background is not varied significantly over  a narrow energy interval.
This method is usually motivated by the fact that the background can be quite complicated and for the observations with large exposure having small statistical errors it may be difficult to find a good background fit, as the data may be e.g. dominated by systematic. While by its nature the sliding window approach can be conservative for a detection of a line (it is more likely to miss a real line  than to detect a fake line), due to  the same reason it is not conservative for the bounds. 
As we try to illustrate below, a formal statistical upper bound obtained by this method may be too strong and exclude a line that can be detected with more precise background modeling. 
If residuals caused by systematic errors are local and absent in most of the energy bins, it may be possible to determine the parameters of the continuum background precisely despite their presence.
A line detected against such a background may still be related to systematic errors, even if it is formally statistically significant.
However, it also can not be excluded as a physical signal. Its nature should therefore be studied by further observations of different targets. 



To demonstrate that the choice of the energy interval used for fitting of the background model can strongly affects the conclusion about the line presence/strength we took the 180'-600' region and modeled the combined
MOS and PN spectra at 3.3-3.8 keV by a sum of two powerlaws: unfolded (that represents the instrumental background) and folded (representing astrophysical emission). For the unfolded component we kept both the index and the normalization at best-fit values found from fitting within our fiducial energy ranges), while for the cosmic emission the powerlaw index and the normalization were allowed to vary. 
On top of these models, the narrow Gaussian line was added. The resulting fit quality was good: χ^2 = 20.7 for 17 d.o.f., null-hypothesis probability is 23.9%. The best-fit line normalization becomes zero, with upper bound being 0.028 cts/cm^2/s/sr (95% CL).  
This upper bound[Remarkably, the bound that we have obtained  from the narrow interval is roughly consistent with the results of <cit.>, even 
though we have used very different (much more primitive) statistical treatment. Indeed, in this exercise we used 
 about 1 Msec of (MOS + PN) observations, while the total exposure of <cit.> is ∼ 30 Msec, Therefore we can roughly re-scale our bound  as ∼ 0.028/√(30)≈[0.005]cts/cm^2/s/sr – close to the value that we have  deduced from the results of <cit.>, see (<ref>).] is stronger than the best-fit line normalization (0.06^+0.02_-0.01 cts/sec/cm^2/sr) reported in Table <ref>, i.e., by doing the sliding window analysis we have excluded the line that we had previously detected with more than 3σ statistical significance.

The reason for that lies in the violation of the basic  assumption of the sliding window method – the background spectrum is not monotonic in our case within statistical errors.
Indeed, as one can see in Table <ref>,  lines at 3.3 and 3.7 keV are present in the spectrum.
By not modeling these lines, we have artificially raised the level of best-fit continuum. 
As the flux in the 3.3 keV line is actually higher than in the 3.5 keV and 3.7 keV lines, this few per cent raise of the continuum level is sufficient to fully absorb two latter lines, strongly decrease the flux in 3.3 keV line (see Fig. <ref>, left panel) and create a stringent upper bound on the flux of 3.5 keV. 
The resulting continuum model, however, strongly over-predicts outside the energy range used for fitting, as Fig. <ref>, right panel, clearly demonstrates.[A similar problem appears in the analysis of <cit.> where two unmodeled lines can be clearly seen on the exclusion plot of Fig. 2 <cit.>.] 

Mathematical meaning of both methods presented in Fig. <ref> is clear, it is up to a physicist to choose which interpretation of the data to use to derive conclusions.

Acknowledgement. We thank B. Safdi for comments on the manuscript. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (GA 694896). DI and OR acknowledge support from the Carlsberg Foundation. The work of DS was supported by the grant for young scientist’s research laboratories of the National Academy of Sciences of Ukraine.

JHEP




    Supplementary Material















§ DATA ANALYSIS AND MODELING
 


We downloaded all publicly available  observations from inner 35^∘ around Sgr A*. We process the observation data using the standard Extended Sources Analysis Software (ESAS), a part of Scientific Analysis Software (XMMSAS v.16.0.0), with an appropriate version of calibration files. The data were cleaned from variable proton component using standard ESAS procedures  and  with default parameters. 
Point sources were automatically detected with the standard SAS procedure . The extraction radius were determined by the distance where surface brightness of the point source equals 0.1 of surrounding background.
We produced source spectra and generated response matrices inside the region overlapping the corresponding annulus using ESAS tasks  and  and the background spectra generated by  and .
For every obtained spectra we calculated the count rate normalized by the value of  keyword. The obtained values were compared inside each of the annulus regions and the high count rate outliers were removed from analysis. 
For the leftover expositions the lightcurve and count rate histogram plots, produced by  and  were visually inspected to ensure the correctness of flared interval automatic removal procedure. Exposures with potential residual soft proton contamination were dropped from further analysis.

For the list of observations in each of the regions Reg1–Reg5 and their cleaned exposures and fields-of-view, see <cit.>.

Individual observation spectra and response files were combined and binned by 45 eV per energy bin (similar to <cit.>) in each of the regions of interest by using the  procedure, a part of  software v.6.24.












Finally, we model the obtained spectra in  spectral package, v. 12.10.0c, also a part of  v.6.24, by using a sum of absorbed cosmic continuum (consisting of thermal and non-thermal component), an instrumental background, and several narrow emission lines of our interest — in addition to the 3.5 keV line, these are the astrophysical emission lines at 3.12, 3.32, 3.68, 3.90 and 4.11 keV. To increase the statistical significance of the lines, we modeled together  spectra from MOS and PN cameras. Similar to our previous works, we used the wide energy range (namely, 2.85-6.0 keV) to model the spectra from regions 1-5, and also included high energies (around E ∼ 10 keV) to better fix the parameters of the instrumental background continuum.









