DETECTION OF AN UNIDENTIFIED EMISSION LINE IN THE STACKED X-RAY SPECTRUM OF GALAXY CLUSTERS

We searched the XMM-Newton archive for galaxy cluster observations that yielded large numbers of X-ray counts. We first selected clusters below a redshift of 0.4; higher-redshift clusters are too faint to contribute significantly to the stacked spectrum. We then calculated the total X-ray counts expected from these XMM-Newton observations using the ROSAT count rates reported in eBCS (Ebeling et al. 2000), NORAS (Böhringer et al. 2000), REFLEX (Böhringer et al. 2004), XBACs (Ebeling et al. 1996), and MACS catalogs (Ebeling et al. 2001) and XMM-Newton exposures. To prevent nearby clusters from dominating the stacked spectrum, we used different cluster count limits for different redshift ranges. We chose clusters with a minimum of 10⁵ counts per cluster for clusters with z < 0.1, and 10⁴ counts per cluster for clusters with redshifts 0.1 < z < 0.4, to have a wide enough range for the redshift-smearing effect. Offset pointings were excluded from the sample. In the end, a sample of 73 clusters was selected. Included in Table 1 are the XMM-Newton observation identification (ObsID) numbers, total MOS and PN clean exposure times, count rates, and our best-fit redshifts (see Section 2.2). The redshift histogram of the sample is given in Figure 1. The count rates reported in Table 1 have been used only for sample selection.

Zoom In Zoom Out Reset image size

The EPIC data processing and background modeling were carried out with the XMM-Newton Extended Source Analysis Software (XMM-ESAS; Kuntz & Snowden 2008; Snowden et al. 2008). We reduced MOS and PN data with the XMM-Newton Science Analysis System (SAS) version 12. Our XMM-Newton analysis is described fully in Bulbul et al. (2012a, 2012b); here we give relevant details.

The light curve filtering was applied to eliminate periods of elevated background. Cleaned events files were created using the good time interval file produced by this process. The net exposure time after filtering the event files for good time intervals is given in Table 1.

Images were created in the 0.4–7.0 keV band for MOS and PN observations and used for detecting point sources with the CIAO tool wavdetect. The images were examined carefully for any missed point sources, as well as for individual MOS CCDs operating in an anomalous state. The CCDs in an anomalous state and all point sources were excluded from further analysis.

Subtle errors in the detector energy gain may cause a fraction of a percent shifts of the location of the emission lines in different X-ray observations of the same cluster. In addition, a redshift measured from the optical observations of a cluster may differ from an X-ray redshift of the gas. To be able to stack spectra in the same frame, we determined the best-fit X-ray redshift for each XMM-Newton observation using the bright Fe lines. These redshifts (Table 1), which correct for both of the above-mentioned effects, were then used to scale the individual spectra in the source frame of each observation, as will be described in Section 2.3. Our selected observations provide adequate statistics to fit an X-ray redshift for each spectrum.

For most clusters, the spectra were extracted within the overdensity radius R₅₀₀. The overdensity radii were calculated using the Vikhlinin et al. (2009) mass–temperature scaling relation for each cluster. Due to the large solid angle of nearby clusters, e.g., Coma, Perseus, and Centaurus, their spectra were extracted within the full field of view (FOV). Redistribution matrix files (RMFs) and ancillary response files (ARFs) were created with the SAS tools rmfgen and arfgen, respectively.

Although we stack the cluster spectra in this work (and end up using only the 2–10 keV band for the line search), it is still important to accurately subtract the background from each individual observation. For each extracted spectrum, we model a superposition of four main background components: quiescent particle background (QPB), cosmic X-ray background emission (including Galactic halo, local hot bubble (LHB), and unresolved extragalactic sources), solar wind CX, as well as residual contamination from soft protons. We use the ROSAT All-Sky Survey (RASS) background spectrum to model the soft X-ray background using the background tool at the High Energy Astrophysics Science Archive Research Center Web site. The RASS spectrum was extracted from an annulus from 1° to 2° surrounding the cluster center, with the assumption that this spectrum reasonably represents the soft X-ray background in the direction of the cluster.

We simultaneously modeled the soft X-ray emission from the LHB or heliosphere with a cool unabsorbed single-temperature thermal component (E ∼ 0.1 keV), while the Galactic hotter halo and intergalactic medium were modeled with an absorbed thermal component (E ∼ 0.2 keV). The energies of the apec model were restricted but allowed to vary with free normalizations. The abundances were set to 1 A_☉. We model the contamination due to unresolved point sources using an absorbed power-law component with a spectral index of α ≃ 1.46 and normalization of 8.88 × 10⁻⁷ photons keV⁻¹ cm⁻² s⁻¹ at ∼1 keV (Kuntz & Snowden 2008). Soft-proton flares are largely removed by the light curve filtering. However, after the filtering some soft-proton residuals may remain in the data and were modeled by including an extra power-law model component and diagonal response matrices provided in the SAS distribution in the final spectral analysis (Snowden et al. 2008).

The EPIC–MOS QPB spectra have two bright instrumental fluorescent lines: the Al–K (1.49 keV) and the Si–K (1.74 keV) lines. The PN QPB spectra have fluorescent lines of Al–K (1.49 keV), Ni–K (7.48 keV), Cu–K (8.05, 8.91 keV), and Zn–K (8.64, 9.57 keV). Since small variations in the gain and the line strengths between the source and background spectra can lead to residuals in the spectral fitting (Kuntz & Snowden 2008) and XMM-ESAS software does not include these instrumental lines in the QPB spectra, we modeled these instrumental lines spectrally by adding Gaussian models to our spectral fits to determine the best-fit energies, widths, and normalizations. The total background was constructed by adding the models for the Al–K, Si–K, Ni–K, Cu–K, and Zn–K lines with the best-fit energies, widths, and normalizations to the QPB produced in the XMM-ESAS analysis for all pointings. These total QBP spectra were directly subtracted from the summed observation to obtain source spectra.

The fitting of the source spectra was done with the spectral fitting package XSPEC 12.8.0 (Arnaud 1996). The 0.3–10 keV energy interval was used for MOS spectra, whereas the 0.4–10.0 keV band was used for the PN fits. To determine the best-fit cluster redshifts for each observation (given in Table 1), the cluster spectra were fit with a standard absorbed multi-temperature collisional equilibrium plasma model (apec; Smith et al. 2001) and AtomDB v2.0.2 (Foster et al. 2012). We did not observe any differences beyond a fraction of a percent in terms of the detector gain variations.

The best way of distinguishing a real spectral feature in a class of distant objects from instrumental artifacts and the X-ray background features is to detect that feature in multiple objects at different redshifts in their rest frame, in which case the line coming from an object will stay at the same energy, unlike the detector artifacts. To accomplish this, we stacked the spectra of our selected 73 clusters, blueshifting them to the source frame using the best-fit X-ray redshift of each observation determined above.

Technically, the energies of the source and background X-ray events were rescaled to the source frame using the best-fit redshifts. The scaled event files were then used to extract the source within r = R₅₀₀ or the full FOV of MOS, and the same extraction region was used for PN observations for nearby clusters that fill the FOV. The particle background spectra were extracted using the scaled filter wheel closed data. Counts from each individual spectrum were co-added into a single stacked spectrum using the FTOOL mathpha to produce the stacked source and the particle background spectra. At the end of the stacking process, we obtained spectra with ∼6 Ms of good cluster exposure with MOS 1 and MOS 2 (that were co-added) and ∼2 Ms with PN for the full XMM-Newton sample.

The RMF and ARF to be used with the stacked spectrum were constructed by averaging the responses for individual observations with proper weighting. The individual RMFs and ARFs were first remapped to the source frame using the best-fit redshifts. The weighing factors for stacking RMFs and ARFs were calculated using the total counts in the energy band we will use for our line search (2–10 keV). These factors (ω_cnt) are given in Table 4. The weighted and normalized ARFs and RMFs were stacked using the FTOOLS addarf and addrmf. These X-ray count-weighted response files were used to model the continuum and the known plasma emission lines; we will also try a different weighting of responses for the possibly nonthermal new line, as will be described below.

For a check, each background-subtracted, blueshifted, single-cluster spectrum was fit with an apec model using the corresponding scaled ARF and RMF to verify that the best-fit redshifts were consistent with zero. For illustration, the co-added MOS and PN source and background spectra of the Perseus Cluster in its source frame are shown in Figure 2. We note that the Fe xxv line is located at its rest energy ∼6.7 keV, while the background and instrumental lines are blueshifted.

Zoom In Zoom Out Reset image size

The stacked MOS and PN source and background spectra of the clusters in the sample are shown in Figure 3. The background spectra show the smearing effect on the background lines, e.g., Al–K (1.48 keV), Si–K (1.75 keV), Cr (5.4 keV), Mn (5.8 keV), Fe–K (6.4 keV), Cu–K (8.05 keV, 8.91 keV), Zn–K (8.64 keV, 9.61 keV), and Au (9.1 keV). They are much less prominent in the stacked spectrum compared with the single-source spectrum shown in Figure 2. Similarly, any residuals from inaccurate background subtraction are smeared. We will see other advantages of this smearing below.

Zoom In Zoom Out Reset image size

After the stacking process, we obtained a total of 8.5 × 10⁶ source counts in the 6 Ms MOS spectra, while the 2 Ms PN stacked spectra have a total of 5.1 × 10⁶ source counts. The line-free apec model with Gaussian lines produces an acceptable fit to the stacked MOS and PN spectra with χ² values of 564.8 for 566 degrees of freedom (dof) (MOS) and 510.5 for 564 dof (PN). After modeling all the known thermal plasma lines in the stacked spectrum, we examined the residuals in each 1 keV band carefully. We found one significant unidentified residual emission feature at E ≈ 3.55–3.57 keV, which is not associated with any plasma emission lines in the band. Near this line, there are four tabulated weak thermal emission lines of K xviii (1s¹ 2s¹ → 1s²) at a rest energy of 3.47 keV, K xviii (1s¹ 2p¹ → 1s²) at 3.51 keV, a dielectronic recombination (DR) line of Ar xvii at 3.62 keV, Ar xvii (1s¹ 3p¹ → 1s²) at 3.68 keV, and K xix (2p¹ → 1s¹) at 3.72 keV.

In order to separate the excess emission feature from these weak contaminating K and Ar lines, we make conservative estimates of their flux using AtomDB. Ideally, line flux measurements would be based on other lines of the same ions; however, there are no other strong K xviii and K xix lines in the spectrum. Therefore, we use the lines from the relatively clean part of the band, namely, the S xvi (2p¹→ 1s¹), Ca xix (1s¹ 2p¹ → 1s²), and Ca xx (2p¹ → 1s¹) lines at 2.63 keV, 3.90 keV, and 4.11 keV, respectively, to estimate the flux of the 3.47 keV, 3.51 keV, 3.68 keV, and 3.72 keV lines. The best-fit flux measurements of these S xvi, Ca xix, and Ca xx lines are given in Table 2.

We assume that the relative abundances of S, Ca, Ar, and K are proportional to their abundances in the solar photosphere (Anders & Grevesse 1989). While this may not be exactly true, it gives a reasonable starting point (we will relax this assumption below). Then, using AtomDB, we calculated the relative emissivity of the K xviii, K xix, and Ar xvii lines compared to the S xvi, Ca xix, and Ca xx lines based on the equilibrium collisional plasma conditions at the various temperatures of our line-free apec components. In practice, the emissivities of K xviii, K xix, and Ar xvii lines are stronger at the lowest temperatures of each model, so the other components can be ignored. The curves in Figure 4 represent the emissivities of K and Ar lines as a function of plasma temperature for the normalizations of the lowest temperature components measured in our spectra.

Zoom In Zoom Out Reset image size

Having obtained the relative theoretical emissivity of the lines from AtomDB, we estimated the flux as

$$\begin{equation} \Gamma _{{ l}} = \Gamma _{{ r}} \sum \limits ^{{ i}} \rm {Norm}_{{ i}}\varepsilon _{{ l}}({ T}_{{ e}})/ \varepsilon _{{ r}}({ T}_{{ e}}), \end{equation} \tag{ 2 }$$

where subscripts l and r represent the lines of interest (K xviii and Ar xvii) and reference lines (S xvi, Ca xix, and Ca xx), respectively, Γ is the flux in the line, ε(T_e) is the calculated emissivity from AtomDB at the electron temperature T_e, and the sum over i represents the different temperature components listed in Table 2 with their normalizations Norm_i. We use 0.1 and 3 times the maximum values of these fluxes as lower and upper bounds, respectively, for the normalizations of the Gaussian lines in the XSPEC fitting. The lower limits of 0.1 is set to avoid the lines vanishing and posing problems for the minimization routine. The factor three represents a conservative allowance for variation of the relative elemental abundances between S and Ca (the measured lines on which the predictions are based) on one hand and K and Ar on the other. (This factor of three is not included in Table 3.)

Since our detected emission line is only 50 eV away from the Ar xvii DR line at rest energy 3.62 keV, we calculated the emissivity of the Ar xvii DR line in a conservative way, using AtomDB v2.0.2. The He-like argon triplet including four lines (known either as w, x, y, z, or R, I1, I2, and F) was summed, since the components cannot be distinguished at the CCD resolution. The two Ar xvii DR lines at 3.62 keV, known in AtomDB as 10077 → 2 and 10078 → 3, and which are the result of an He-like Ar ion recombining to Li-like Ar and emitting a photon at 3.62 keV, were similarly extracted and summed. The right panel of Figure 4 shows the comparison of the emissivity of Ar xvii DR and the He-like argon triplet at E ≈ 3.12 keV. To model the flux of the Ar xvii DR line in our spectral fits in a conservative way, we set the lower and upper limits of the flux to be 0.001 and 0.01 times the flux of the He-like Ar, respectively. The upper limit corresponds to the highest flux that Ar xvii DR can have for the ICM plasma temperatures that we see in our spectra (this will be further discussed in Section 3.4). For a plasma temperature of ∼2 keV, the lowest temperature observed in our samples, the ratio of the flux of the Ar xvii DR line to that of the He-like Ar line corresponds to 0.01 and thus was chosen as an upper limit (see Figure 4, right panel). The lower limit has been set to avoid problems with the fitting procedure.

Once the lower and upper limits on flux estimates of the K xviii, K xix, and Ar xvii lines were set, we performed the fit in a narrower 3–6 keV energy band (to avoid strong S and Si lines below 3 keV and Fe lines above 6 keV). This band is sufficiently wide to measure the continuum accurately (to better than 1%). The weak residual emission line at E ≈ 3.57 keV was detected in the fits. The excess emission after the Gaussian K and Ar lines were included in the model at their maximum fluxes (as described above) in MOS and PN spectra is shown in Figure 5. Figure 6 shows the excess in the rebinned MOS spectrum of the full sample. We then added a Gaussian model to fit the remaining residuals, leaving its flux and energy to vary. The fit was improved by Δχ² of 22.8 for MOS and Δχ² of 13.9 for PN for an additional 2 dof (energy and normalization). The best-fit energy of the added Gaussian line is 3.57 ± 0.02 (0.03) keV in the stacked MOS and 3.51 ± 0.03 (0.04) keV in the stacked PN observations. The line energies from MOS and PN are in significant tension, 2.8σ apart (Figure 9). However, given the systematic uncertainties of the fitting procedure, we consider it acceptable; this tension disappears once another level of complexity is introduced in modeling (see Section 3.5 below). The width of the new line is unresolved and broadened only by the instrumental response. This is the only significant unidentified feature we have detected in the 2–10 keV band of MOS and PN spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To measure the flux of this line, we have to use a statistically proper response file, which will depend on the physical interpretation of the line. If the line were coming from the thermal plasma, then the same spectral responses that were used for the thermal components are appropriate. However, there are no known thermal plasma lines at this energy, so we explore a possible interpretation of the detected line as a decay signature of the sterile neutrino (see Section 1). In this interpretation, the spectral fitting procedure has to be slightly modified. In particular, when co-adding the instrumental responses used for the dark matter (DM) line component, the individual cluster responses should be weighted by the factor ω_dm proportional to the estimated dark matter photon flux from each cluster (as opposed to the X-ray flux used for the response averaging so far).

These response files will be solely used to measure the flux of the detected 3.57 keV line; for the rest of the components, clearly originating in the ICM, the X-ray flux weighting is correct. The dark matter response weighting was done using the following approach.

The surface brightness of the DM decay signal is proportional to the DM column density S_DM = ∫_losρ_DM(r)dr. The observed photon flux from the DM decay into a solid angle Ω_FOV is given by

$$\begin{equation} F_{{\rm DM}} = \frac{M_{{\rm DM}}^{{\rm FOV}} }{4\pi D_{L}^{2} } \frac{\Gamma _{\gamma }}{m_{s}} (1+z)\ \ \ \rm {photons\ cm^{-2}\ s^{-1}}, \end{equation} \tag{ 3 }$$

where Γ_γ and m_s are the decay rate and mass of the sterile neutrino (see Equation (1) and Pal & Wolfenstein (1982)), $M_{{\rm DM}}^{{\rm FOV}}$ is the projected DM mass within the spectral extraction region (R_ext, which is either R₅₀₀ or R_FOV), and D_L is the luminosity distance.

The DM mass projected along the line of sight is

$$\begin{equation} M_{{\rm DM}}^{{\rm FOV}} = \int _{{\rm los}} \rho _{{\rm DM}} (r) dr, \end{equation} \tag{ 4 }$$

where ρ_DM(r) is the distribution of dark matter determined by the Navarro–Frenk–White (NFW) profile (Navarro et al. 1997) and given by

$$\begin{equation} \rho _{{\rm DM}}(r) = \frac{\rho _{c}}{(r/r_{s})(1+r/r_{s})^{2}}, \end{equation} \tag{ 5 }$$

where ρ_c is a characteristic density and r_s is a scale radius. The integration of the dark matter distribution within the extraction radius (given in Table 4) is along the line of sight. An extraction radius of 700'' was used for the clusters larger than the FOV of XMM-Newton, while an extraction radius of R₅₀₀ was used for the clusters smaller than the FOV.

Table 4. Weighting Factors Used for Stacking Response Matrices

Cluster	$M^{{\rm proj}}_{{\rm DM}}$ (Rext)	Rext	$M^{{\rm proj}}_{{\rm DM}}$/D2	Exp/Exptot	ωdm	ωcnt
	(1014 M☉)	(Mpc)	(1010 M☉ Mpc−2)
	(1)	(2)	(3)	(4)	(5)	(6)
Centaurus	0.63	0.17	2.41	0.044	0.139	0.074
A1060	0.59	0.21	1.82	0.010	0.024	0.009
A262	0.52	0.24	1.24	0.015	0.025	0.011
Perseus	1.49	0.24	2.89	0.048	0.181	0.39
AWM7	0.86	0.24	1.82	0.045	0.106	0.061
Coma	2.72	0.33	2.78	0.026	0.094	0.06 2
A3581	1.32	0.27	1.35	0.028	0.050	0.013
Ophiuchus	4.14	0.38	3.05	0.009	0.037	0.032
A4038	1.31	0.39	0.91	0.008	0.010	0.007
A496	2.29	0.44	1.24	0.038	0.061	0.044
A2063	1.92	0.48	0.88	0.008	0.009	0.0057
A2147	2.06	0.47	0.96	0.002	0.003	0.0016
A3571	3.94	0.53	1.42	0.007	0.013	0.012
A3558	3.40	0.64	0.82	0.012	0.013	0.01
A4059	2.75	0.61	0.27	0.010	0.003	0.007
Triangulum Australis	7.58	0.68	1.66	0.002	0.006	0.003
Hydra A	2.68	0.72	0.51	0.025	0.016	0.023
A754	7.91	0.72	1.48	0.004	0.008	0.0032
A2319	6.93	0.72	1.31	0.003	0.0004	0.033
Cygnus A	3.81	0.72	0.72	0.005	0.005	0.004
AS1101	1.95	0.79	0.34	0.025	0.011	0.0136
A3112	4.44	0.96	0.45	0.054	0.034	0.0337
A2597	3.61	0.96	0.29	0.004	0.002	0.002
A478	8.30	1.10	0.61	0.019	0.014	0.017
PKS0745−19	10.03	1.27	0.52	0.005	0.003	0.003
A2811	5.29	1.08	0.15	0.007	0.001	0.0018
A2034	8.07	1.25	0.35	0.005	0.002	0.002
RXC J0616.8−4748	3.97	0.99	0.16	0.0069	0.001	0.0007
RXC J0145.0−5300	6.11	1.14	0.25	0.011	0.003	0.0025
RXC J1044.5−0704	3.05	0.89	0.09	0.007	0.0009	0.0014
A1068	4.44	1.01	0.12	0.005	0.0009	0.0012
RXC J2218.6−3853	6.68	1.16	0.20	0.005	0.001	0.0013
RXC J0605.8−3518	4.91	1.05	0.14	0.005	0.001	0.0013
A1413	9.09	1.29	0.24	0.053	0.016	0.018
A2204	8.86	1.27	0.21	0.019	0.005	0.010
A3888	8.57	1.26	0.20	0.013	0.003	0.004
RXC J0958.3−1103	6.58	1.15	0.15	0.002	0.0005	0.0006
A545	10.79	1.36	0.25	0.002	0.0005	0.0004
RXC J2014.8−2430	6.18	1.12	0.13	0.006	0.001	0.002
RX J1720.1+2638	6.64	1.25	0.13	0.016	0.003	0.004
RXC J0645.4−5413	8.55	1.47	0.16	0.005	0.001	0.001
A1201	5.78	1.10	0.11	0.015	0.002	0.0017
A1914	13.93	1.67	0.26	0.006	0.002	0.002
A2345	7.65	1.20	0.14	0.015	0.003	0.002
A2218	7.48	1.19	0.13	0.012	0.002	0.002
A2254	7.47	1.19	0.13	0.017	0.003	0.002
A665	9.50	1.29	0.16	0.006	0.001	0.0015
A1689	12.55	1.42	0.20	0.010	0.002	0.004
A383	3.48	0.92	0.05	0.008	0.0005	0.009
A520	7.75	1.20	0.11	0.009	0.001	0.001
A2163	26.34	1.80	0.34	0.003	0.001	0.001
A209	8.82	1.25	0.11	0.005	0.0007	0.0007
A963	6.81	1.15	0.07	0.006	0.0006	0.001
RXC J1504.1−0248	8.87	1.25	0.09	0.011	0.003	0.004
MS 0735.6+7421	3.89	0.95	0.04	0.014	0.0008	0.001
A773	9.34	1.27	0.11	0.004	0.0005	0.0004
AS0592	13.27	1.42	0.14	0.008	0.002	0.0017
A2390	12.07	1.38	0.13	0.003	0.0005	0.0008
A2667	9.66	1.28	0.10	0.006	0.0007	0.0011
A267	4.83	1.01	0.05	0.002	0.0001	0.0005
RXC J2129.6+0005	3.06	0.87	0.03	0.0097	0.0004	0.001
RXC J1314.4−2515	8.61	1.22	0.07	0.010	0.0009	0.004
A1835	12.15	1.37	0.10	0.037	0.005	0.009
A1758	4.54	1.04	0.03	0.009	0.0004	0.0008
A1763	10.47	1.32	0.11	0.004	0.0005	0.0005
A689	22.51	1.66	0.15	0.002	0.0001	0.0001
ZW 3146	6.72	1.11	0.04	0.059	0.003	0.010
A781	5.57	1.04	0.03	0.018	0.0007	0.001
Bullet	15.24	1.45	0.09	0.006	0.0007	0.001
MS 2137.3−2353	4.31	0.95	0.02	0.003	0.0001	0.0002
MACS J2229.7−2755	3.51	0.88	0.02	0.009	0.0001	0.0006
MACS J1532.8+3021	4.85	0.97	0.02	0.003	0.0007	0.0003
AS1063	16.80	1.48	0.07	0.004	0.0004	0.0008

Notes. Columns 1 and 2 show the estimated projected dark matter masses in the spectral extraction radii $M_{{\rm DM}}^{{\rm proj}}$ (R_ext) and the extraction radii R_ext in Mpc; Column 3 is the projected dark matter masses per distance squared; Column 4 shows the ratio of the exposure time to the total exposure stacked for each cluster; Column 5 is the weighting factors (ω_dm) calculated based on the predicted dark matter flux used in the stacking of ARFs and RMFs of each cluster in the sample. These stacked ARFs and RMFs were then used to determine the flux of the detected line, and Column 6 shows the weighting factors (ω_cnt) calculated based on the total counts in the fitting band. The response files that were stacked using these factors were utilized to model plasma emission lines.

Download table as:  ASCIITypeset images: 1 2

The expected contribution of each cluster i to the total DM line flux in the stacked spectrum is

$$\begin{equation} \omega _{i,{\rm dm}} = \frac{M_{i, {\rm DM}}^{{\rm proj}}(<R_{{\rm ext}})(1+z_{i})}{4\pi D_{i,L}^{2}}\ \frac{e_{i}}{e_{{\rm tot}}}, \end{equation} \tag{ 6 }$$

where z_i is the redshift of the ith cluster and e_i and e_tot are the exposure time of the ith cluster and the total exposure time of the sample, respectively.

The dark matter mass within the extraction radius is estimated as

$$\begin{equation} M_{{\rm DM}}(R_{{\rm ext}}) = M_{{\rm tot}} (R_{{\rm ext}}) - M_{{\rm gas}}(R_{{\rm ext}})-M_{*}(R_{{\rm ext}}), \end{equation} \tag{ 7 }$$

where M_tot(R_ext), M_gas(R_ext), and M_*(R_ext) are the total mass, gas mass, and stellar mass in the extraction radius R_ext, respectively. The observed Vikhlinin et al. (2009) temperature–mass scaling relation was used to infer total masses for the intracluster gas temperatures measured from the XMM-Newton observations. The gas mass is determined following the method described in Bulbul et al. (2010). The contribution of stars to the total baryon budget is modest at large radii but more important in the cluster centers because of the presence of cD galaxies. At large radii (⩾R₅₀₀), M_* is 10%–15% of the gas mass (Lin & Mohr 2004; Vikhlinin et al. 2006). Stellar masses of each cluster were determined using the stellar mass–total mass scaling relation (Gonzalez et al. 2013). The calculated dark matter masses were corrected using this factor. The projected dark matter masses within R_ext were then determined by projecting NFW profiles (Bartelmann 1996; Golse & Kneib 2002; Loewenstein et al. 2009). We used a concentration parameter c₅₀₀ = 3 from the Vikhlinin et al. (2006) c − M₅₀₀ scaling relation and the median total mass within R₅₀₀ of the full sample, which is ∼6 × 10¹⁴ M_☉. The projected dark matter mass within each spectral extraction radius is given in Table 4.

Weights for the responses to be included in the stacked-spectrum response were calculated as follows. The number of dark matter decay photons in each cluster spectrum is

$$\begin{equation} S_{i} = \alpha \ \omega _{i,{\rm dm}}\ e_{{\rm tot}} \ A_{i}, \end{equation} \tag{ 8 }$$

where A_i is the ancillary response (the instrument effective area) at photon energy E/(1 + z_i), and α is the ratio of the decay rate of sterile neutrinos to the sterile neutrino mass m_s (here we denote α ≡ Γ_γ/m_s). The total number of dark matter photons in the stacked line is

$$\begin{eqnarray} S_{{\rm line}}&=&\sum _{i=0}^{i=73}S_{i} \nonumber\\ &=& \alpha \ \omega _{{\rm tot}} \ e_{{\rm tot}}\ A_{\omega }, \end{eqnarray} \tag{ 9 }$$

where the weighted ARF A_ω is a function of the total weight ω_tot,

$$\begin{equation} A_{\omega }= \sum _{i} \frac{\omega _{i}}{\omega _{{\rm tot}}}A_{i}, \end{equation} \tag{ 10 }$$

and

$$\begin{equation} \omega _{{\rm tot}}= \sum _{i}\omega _{i}. \end{equation} \tag{ 11 }$$

The weighted responses A_ω were used to model our new line, while X-ray count-weighted response files were used to model the other known emission lines and the continuum components.

For MOS, the flux in the 3.57 keV line was $4.0^{+0.8}_{-0.8}$ ($^{+1.8}_{-1.2}$) $\times 10^{-6} \ \rm {photons}\ \rm {cm}^{-2}\ \rm s^{-1}$, where the errors are 68% (90%). For PN, at the best-fit energy of 3.51 keV, the line flux is $3.9^{+0.6}_{-1.0}$ ($^{+1.0}_{-1.6})\,\times 10^{-6} \ \rm {photons}\ \rm {cm}^{-2}\ \rm s^{-1}$. If we fix the line energy from the MOS fit, for PN we obtain the flux $2.5^{+0.6}_{-0.7}$ ($^{+1.0}_{-1.1})\,\times 10^{-6} \ \rm {photons}\ \rm {cm}^{-2}\ \rm s^{-1}$.

We note that the line energy detected in the stacked PN observations of the full sample is consistent with the K xviii line at 3.515 keV. However, the measured flux from this line is a factor of 20 above the expected flux of the K xviii line, estimated based on the measured fluxes of the S xvi, Ca xix, and Ca xx lines and assuming a consistent relative abundance for K xviii along with the plasma temperature from AtomDB. In addition, the detected energy in the stacked MOS observations of the full sample is 3.5σ away from the K xviii line. This will be further discussed later.

Since this is a blind search, in which the line energy is not known a priori, to estimate the significance of our detection, we must take into account the "look elsewhere" effect. We have examined ∼70 independent energy resolution elements in our search band and would accept a line detection in any of these bins. Taking this into account, our 4σ–5σ detections correspond to the probability of falsely detecting a line at an unknown energy of 0.004% for MOS and 0.4% for PN. However, the line is found at a consistent energy (or at least in the same independent resolution element) in these two completely independent samples coming from different instruments. The statistical chance of such a false detection at the same energy is negligibly low. We caution that these are just the rough estimates of the statistical probabilities; systematic uncertainties are also important (Section 6).

Because estimating statistical significance of faint line features is a notoriously ill-behaved problem, we have verified the above estimate with Monte Carlo simulations. We used the PN detection for this test, because its significance is lower and a Monte Carlo estimate can be done using a reasonable number of trials. We generated 1000 random realizations of a spectrum using a model in XSPEC with no extra emission line and fit each of them with a model that included an additional line at an arbitrary location and flux. We then counted the realizations in which the model with the additional line improved the fit by Δχ > 11.2, which corresponds to our PN detection. This false detection occurred in 4 cases out of 1000, in agreement with the above 0.4% probability of false detection in the stacked PN spectrum.

We also fit the same MOS and PN spectra using the X-ray count-weighted responses, to check if the detection is dependent on the response weighting. For MOS, the flux of the detected line was $4.1^{+1.0}_{-0.9}$ ($^{+1.8}_{-1.6}$) $\times 10^{-6} \ \rm {photons}\ \rm {cm}^{-2}\ \rm s^{-1}$; the fit was improved by Δχ² of 21.8 for 2 dof. For PN, the line flux was $3.9^{+1.3}_{-1.0}$ ($^{+2.1}_{-2.0})\,\times 10^{-6} \ \rm {photons}\ \rm {cm}^{-2}\ \rm s^{-1}$, while the fit was improved by Δχ² of 13.8 for 2 dof. This shows that the detection is robust and the flux is independent of the response scaling.

We will discuss the possible physical interpretations of this emission line in Section 5. Here we will push forward with one possible interpretation of the detected line, sterile neutrino decay, because we need to describe the calculation of certain quantities that will be used below for cross-checks and comparison of the subsamples of our full sample.

For a DM particle decaying radiatively with E_γ = m_s/2, the detected flux from a clump of matter of a known mass can be converted into the decay rate. The energy of the detected line corresponds to a sterile neutrino particle mass of m_s = 7.1 ± 0.07 keV, assuming that the dark matter is solely composed of sterile neutrinos. The relation between the flux and mass implies a mixing angle of

$$\begin{eqnarray} {\rm sin}^{2} (2\theta) &=& \frac{F_{{\rm DM}}}{12.76\ \rm {cm^{-2}} \ s^{-1}} \left(\frac{\rm 10^{14}\ { M}_{\odot }}{M_{{\rm DM}}^{{\rm FOV}}}\right) \nonumber\\ &&\times \left(\frac{D_{L}}{100\ \rm {Mpc}}\right)^{2} \left(\frac{1}{1+z}\right) \left(\frac{1\ \rm {keV}}{m_{s}}\right)^{4}, \end{eqnarray} \tag{ 12 }$$

where F_DM is the observed DM flux.

Using the ω_dm and the projected dark matter masses given in Table 4, we find that the weighted projected dark matter mass per distance squared is 1.82 × 10¹⁰ M_☉ Mpc⁻² for the full sample observed with XMM-Newton MOS. Using Equation (3), one can calculate the mixing angle for the full MOS cluster sample to be $\sin ^{2}(2\theta)= 6.8^{+1.4}_{-1.4}\ (^{+2.0}_{-3.0})\ \times 10^{-11}$. The PN observations of the full sample give a mixing angle measurement of $\sin ^{2}(2\theta)= 6.7^{+1.7}_{-1.0}\ (^{+2.7}_{-1.7})\ \times 10^{-11}$ for a weighted mass per distance squared of 1.80 × 10¹⁰ M_☉ Mpc⁻². These are given in Table 5. The PN and MOS full-sample measurements are consistent with each other and the constraints placed by previous studies, e.g., the unresolved cosmic X-ray background (CXB) in the Chandra Deep Fields (Abazajian et al. 2007) and the XMM-Newton blank-sky background spectrum (Boyarsky et al. 2006), Chandra observations of the Milky Way (Riemer-Sørensen et al. 2006), Chandra observation of the Bullet Cluster (Boyarsky et al. 2008), Chandra observations of the dwarf galaxy Draco (Riemer-Sørensen & Hansen 2009), and XMM-Newton limits from M31 and Willman 1 and Fornax dwarf galaxies (Boyarsky et al. 2010; Watson et al. 2012), as shown in Figure 13(a). It is in marginal (∼90% significance) tension with the most recent Chandra limit from M31 (Horiuchi et al. 2014), as shown in Figure 13(b).

Table 5. Properties of the Unidentified Emission Line

Sample		Inst.	Energy	Flux	χ2	Δχ2	$M_{DM}^{{\rm proj}}$/D2	sin2(2θ)
			(keV)	(10−6 photons cm−2 s−1)	(dof)	(Δ dof)	(1010 M☉ Mpc−2)	(10−11)
		(1)	(2)	(3)	(4)	(5)	(6)	(7)
		MOS	3.57 ± 0.02 (0.03)	4.0 $^{+0.8}_{-0.8}$ $(^{+1.8}_{-1.2}$)	564.8	22.8	1.82	6.8 $^{+1.4}_{-1.4}$ ($^{+2.0}_{-3.0}$)
					(566)	(2)
Full	XMM
sample		PN	3.51 ± 0.03 (0.05)	3.9 $^{+0.6}_{-1.0}$ ($^{+1.0}_{-1.6}$)	510.5	13.9	1.80	6.7 $^{+1.7}_{-1.0}$ ($^{+2.7}_{-1.7}$)
					(564)	(2)
		PN	3.57⋆	2.5 $^{+0.6}_{-0.7}$ ($^{+1.0}_{-1.1}$)	510.5	11.2	1.80	4.3 $^{+1.2}_{-1.0}$ ($^{+1.8}_{-1.7}$)
					(564)	(1)
		MOS	3.57⋆	15.9 $^{+3.4}_{-3.8}$ ($^{+6.7}_{-5.5}$)	562.3	17.1	2.68	18.2 $^{+4.4}_{-3.9}$ ($^{+12.6}_{-11.5}$)
Coma +					(569)	(1)
Centaurus +	XMM
Ophiuchus		PN	3.57⋆	<9.5	377.8	...	...	<10.9
					(387)
		MOS	3.57⋆	21.4 $^{+7.0}_{-6.3}$ ($^{+11.2}_{-10.5}$)	596.1	12.8	2.82	23.3 $^{+7.6}_{-6.9}$ ($^{+12.2}_{-11.5}$)
Perseus					(574)	(1)
(without	XMM
the core)		PN	3.57⋆	<16.1	539.1	...	...	<17.6
					(553)
		MOS	3.57⋆	52.0 $^{+24.1}_{-15.2}$ ($^{+37.0}_{-21.3}$)	613.8	15.7	2.89	55.3 $^{+25.5}_{-15.9}$ ($^{+39.3}_{-22.6}$)
Perseus					(574)	(1)
(with	XMM
the core)		PN	3.57⋆	<17.7	539.4	...	...	<18.8
					(554)
		MOS	3.57⋆	2.1 $^{+0.4}_{-0.5}$ ($^{+0.8}_{-0.8}$)	547.2	16.5	1.08	6.0 $^{+1.1}_{-1.4}$ ($^{+2.3}_{-2.3}$)
All					(573)	(1)
other	XMM
clusters		PN	3.57⋆	2.0 $^{+0.3}_{-0.5}$ ($^{+0.5}_{-0.8}$)	741.9	15.8	1.15	5.4 $^{+0.8}_{-1.3}$ ($^{+1.3}_{-2.1}$)
					(751)	(1)
		ACIS-S	3.56 ± 0.02 (0.03)	10.2 $^{+3.7}_{-3.5}$ ($^{+4.8}_{-4.7}$)	201	11.8	0.72	40.1 $^{+14.5}_{-13.7}$ ($^{+18.9}_{-18.2}$)
					(197)	(2)
Perseus	Chandra
		ACIS-I	3.56⋆	18.6 $^{+7.8}_{-8.0}$ ($^{+12.0}_{-16.0}$)	152.6	6.2	1.86	28.3 $^{+11.8}_{-12.1}$ ($^{+18.2}_{-24.3}$)
					(151)	(1)
Virgo	Chandra	ACIS-I	3.56⋆	<9.1	189.1	...	2.41	<10.5
					(155)

Notes. Columns 2 and 3 are the measured rest energy and flux of the unidentified line in the units of photons cm⁻² s⁻¹ at the 68% (90%) confidence level. The energies with asterisks are frozen to the indicated values; Columns 4 and 5 show the χ² before the line is added to the total model and change in the χ² when an additional Gaussian component is added to the fit; Column 6 is the weighted ratio of mass to distance squared of the samples; and Column 7 shows the mixing angle limits measured in each sample. Reported constraining limits are 90% confidence upper limits. Energies marked with star symbols were held fixed during the model fitting.

Download table as:  ASCIITypeset image

For the PN flux for the line fixed at the best-fit MOS energy, the corresponding mixing angle is $\sin ^{2}(2\theta)= 4.3 ^{+1.2}_{-1.0}\ (^{+1.8}_{-1.7})\ \times 10^{-11}$. This measurement is consistent with that obtained from the stacked MOS observations at a 1σ level. Since the most confident measurements are provided by the highest signal-to-noise ratio stacked MOS observations of the full sample, we will use the flux at energy 3.57 keV when comparing the mixing angle measurements for the sterile neutrino interpretation of this line.

We now divide the full cluster sample into three independent subsamples in order to check that our line does not originate from any single object. The full stacked spectra examined in Section 3.1 have a significant contribution of photons from several nearby bright clusters, e.g., Perseus, Coma, Centaurus, and Ophiuchus. In order to determine whether the line detection is dominated by these bright sources, we excluded them from the sample and stacked the MOS and PN spectra of the remaining 69 fainter galaxy clusters. We have performed the stacking process following the same approach described in Section 2.3. A total of 4.9 Ms of good stacked MOS and 1.7 Ms good stacked PN exposure were obtained for this subsample. The weighted mean redshift was 0.06. The stacked MOS and PN spectra contain 34% (2.95 × 10⁶ source counts) and 55% (2.79 × 10⁶ source counts) of the total source counts of the full cluster sample.

We fit the stacked spectra using the line-free apec model and additional Gaussian models as described in Section 3.1 in the 3–6 keV band. The best-fit temperatures, normalizations of the line-free apec model, and the fluxes of S xvi, Ca xix, and Ca xx lines are given in Table 2. We then carefully examined the spectra for any unidentified emission features in the 3.4–3.7 keV energy interval. Similarly, we determined the maximum fluxes of the K xviii, K xix, and Ar xvii lines based on the plasma temperatures and fluxes of hydrogen-like S xvi, helium-like Ca xix, and hydrogen-like Ca xx lines at 2.63 keV, 3.90 keV, and 4.11 keV, measured from the spectral fits, and AtomDB as described in Section 3.1. As before, the lower and upper limits of the fluxes of K xviii, K xix, and Ar xvii lines were set to 0.1−3 times the maximum predicted fluxes. The Ar xvii DR line flux at 3.62 keV was allowed to vary between 10⁻³ and 10⁻² of the Ar xvii triplet line at 3.12 keV.

We obtained an acceptable fit to the stacked MOS spectrum of these 69 clusters. The total χ² was 557 for 573 dof. Adding in an extra Gaussian model to the MOS spectrum at 3.57 keV improved the fit by Δχ² of 16.5 for an additional degree of freedom. We found that the best-fit flux was 2.1 $^{+0.4}_{-0.5}$ ($^{+0.8}_{-0.8}$) × 10⁻⁶ photons cm⁻² s⁻¹. This flux corresponds to a mixing angle of sin ²(2θ) = 6.0 $^{+1.1}_{-1.4}$ ($^{+2.3}_{-2.3}$) × 10⁻¹¹, consistent with the mixing angle estimates obtained from the full sample.

The overall fit to the stacked PN spectrum for these 69 clusters was acceptable with a total χ² of 741.9 for 751 dof. Adding an extra Gaussian line at 3.57 keV improved the fit by Δχ² of 15.8 for an additional degree of freedom. The PN spectrum yields the best-fit flux detection of 2.0 $^{+0.3}_{-0.5}$ ($^{+0.5}_{-0.8}$) $\times 10^{-6}\ \rm {photons}\ \rm {cm}^{-2}\ \rm s^{-1}$. The mixing angle obtained from the stacked PN observations sin ²(2θ) = 5.4 $^{+0.8}_{-1.3}$ ($^{+1.3}_{-2.1}$) × 10⁻¹¹ is also consistent with the estimates from the full sample. Bottom panels in Figure 5 show the residuals before and after a Gaussian line is added at 3.57 keV to MOS and PN spectral fits.

We now check the MOS and PN spectra of the three dominant nearby clusters, Coma, Ophiuchus, and Centaurus. A total of 525.3 ks of good stacked MOS and 168 ks good stacked PN exposure times were obtained for this subsample. The total source counts obtained in the MOS and PN spectra were 3.2 × 10⁶ and 2.1 × 10⁶, respectively.

We performed the fits as above. The best determinations for the continuum temperature and normalizations and the fluxes of the S xvi, Ca xix, and Ca xx are given in Table 2. We detected an excess emission feature in the same band, i.e., 3.4–3.7 keV as in the stacked MOS spectra. To determine the flux of the emission line at 3.57 keV, we estimated the maximum fluxes of the K xviii, K xix, and Ar xvii lines using the AtomDB and the measured fluxes of S xvi, Ca xix, and Ca xx as described in Section 3.1. Using 0.1 and 3 times these fluxes as lower and upper limits, we found that the unidentified line has a flux of 1.6$^{+0.3}_{-0.4}$ ($^{+0.7}_{-0.6}$) × 10⁻⁵ photons cm⁻² s⁻¹ in the stacked MOS observations. Adding this Gaussian to the model improves the fit by Δχ² of 17.1 for an additional degree of freedom for the stacked MOS spectrum.

We then allowed the energy of the additional Gaussian model to vary to test whether the energy measured from two different samples is the same. The best-fit energy obtained from the stacked MOS observations of Coma, Centaurus, and Ophiuchus clusters was 3.56 ± 0.02 (0.03), with a flux of 1.6$^{+0.52}_{-0.44}$ ($^{+0.81}_{-0.70}$) × 10⁻⁵ photons cm⁻² s⁻¹. This measurement is consistent with the energy measured in the MOS observations of the full sample. The sterile neutrino mixing angle that corresponds to this flux is $\sin ^{2} (2\theta) =18.2^{+4.4}_{-3.9} \ (^{+12.6}_{-11.5})\ \times 10^{-11}$, consistent at 2σ with the full-sample value.

The fits to the stacked PN observations did not need an additional Gaussian line and resulted in a non-detection. This could be due to the low count statistics of the stacked PN observations (168 ks clean time). A 90% upper limit on the flux of this line at 3.57 keV is 9.5 × 10⁻⁶ photons cm⁻² s⁻¹ from this spectrum; the upper limit on the mixing angle from this flux limit is consistent with the full-sample and MOS detections.

Initially, we extracted the spectrum of the Perseus Cluster using the entire MOS and PN FOV. We have co-added the XMM-Newton MOS and PN observations of the Perseus Cluster in the cluster's frame. The total exposure time in the stacked MOS spectrum was 317 ks with a total of 7 × 10⁶ source counts in the 2–10 keV band and 38 ks total exposure with 2 × 10⁶ source counts in the stacked PN observations.

Following the same approach we used for modeling the full cluster sample, we first fit the MOS and PN observations with the line-free apec model and additional Gaussian models. Count-weighted responses were used to fit the plasma emission lines and the continuum emission. Probing the 3–4 keV band, the MOS observations revealed residuals around 3.57 keV, at the same energy band where we detected line emission in the previous samples. The left panel of Figure 7 shows the detection in the co-added MOS observations of the Perseus Cluster. Using the limits on the K and Ar lines (Table 3) as above and adding a Gaussian model to the MOS spectrum at the fixed energy of 3.57 keV improved the fit by Δχ² of 15.7. The best-fit flux at 3.57 keV was 5.2$^{+2.4}_{-1.5}$ ($^{+3.7}_{-2.1}$) × 10⁻⁵ photons cm⁻² s⁻¹.

Zoom In Zoom Out Reset image size

This flux corresponds to a mixing angle of $\sin ^{2} (2\theta)=5.5^{+2.6}_{-1.6} \ (^{+3.9}_{-2.3})\ \times 10^{-10}$. This angle not only is an outlier in our measurements from the other samples but is also not consistent with the upper limits on the mixing angle at this value of m_s from the previous studies (e.g., Horiuchi et al. 2014).

We were unable to detect the line in the short (38 ks clean time) PN observation of Perseus and placed a 90% upper limit on the flux of the line of 17.7 photons cm⁻² s⁻¹, which corresponds to an upper limit of sin ²(2θ) < 1.9 × 10⁻¹⁰, consistent with the MOS detection. Figure 7 shows both XMM-Newton Perseus spectra.

Since this is a single-cluster spectrum, we first check whether the Perseus signal is not an artifact of our blueshifting procedure. For this we fit the original, redshifted MOS spectrum with a line-free apec model. We obtained a best-fit χ² of 463 for 385 dof. Adding a Gaussian line at 3.57 keV (rest energy) improved the fit by Δχ² of 16 for an additional degree of freedom. The best-fit flux was 5.3 ± 1.2 (2.0) × 10⁻⁵ photons cm⁻² s⁻¹, which is in agreement with the flux obtained from the blueshifted spectrum. We conclude that our detection is independent of shifting the spectrum.

Not ready to abandon the sterile neutrino explanation based on the line flux incorrectly scaling with cluster mass that we see for Perseus, we tried to investigate possible astrophysical reasons behind the excess of the line flux in Perseus. First, we investigated the dependence of the energy and flux of this unidentified line on the AtomDB predicted fluxes of nearby lines, i.e., the K xviii line at 3.51 keV and the Ar xvii DR line at 3.62 keV. Allowing the energy of the Gaussian component to vary produced a best fit for an energy of 3.56 $^{+0.01}_{-0.02}$ ($^{+0.02}_{-0.03}$) keV, with a flux of 6.0$^{+1.8}_{-1.4}$ ($^{+2.4}_{-1.7}$) × 10⁻⁵ photons cm⁻² s⁻¹ (χ² of 598.1 for 572 dof). The best-fit energy is consistent with the energy measured from the MOS observations of the full sample. However, the fluxes of the nearby K xviii line at 3.51 keV and the Ar xvii DR line at 3.62 keV were at their allowed upper limits predicted from AtomDB. Relaxing the upper limits has shifted the line energy higher, to 3.59 $^{+0.01}_{-0.03}$ ($^{+0.02}_{-0.04}$) keV, with a flux of 5.5$^{+1.7}_{-0.8}$ ($^{+3.7}_{-1.5}$) × 10⁻⁵ photons cm⁻² s⁻¹ giving a slightly better fit (χ² of 594.5 for 572 dof). We note that the line energy of this extra line gets close to the Ar xvii DR line at 3.62 keV. So we removed the extra Gaussian line and re-fit the Perseus spectrum removing the upper limits on the Ar xvii DR line. We obtained only a slightly worse fit than the previous case, with a χ² of 598.8 (574 dof). The measured flux of the Ar xvii DR line at 3.62 keV in this case was 4.8$^{+0.7}_{-0.8}$ ($^{+1.3}_{-1.4}$) × 10⁻⁵ photons cm⁻² s⁻¹, which is a factor of 30 above the predicted maximum flux of the Ar xvii DR line based on the measured flux of the Ar xvii line at ∼3.12 keV and AtomDB line rates. The predicted maximum flux of the Ar xvii DR line for the Perseus spectrum was 1.6 × 10⁻⁶ photons cm⁻² s⁻¹ (<0.01 times the flux of the Ar xvii triplet at ∼3.12 keV).

This test showed that the line detected in the Perseus Cluster could also be interpreted as an abnormally bright Ar xvii DR line. We note, however, that obtaining such a bright DR line relative to the He-like triplet at 3.12 keV is problematic. The emissivity of the satellite line peaks at kT = 1.8 keV and declines sharply at lower temperatures, in addition to the change in the ionization balance, which reduces the Ar⁺¹⁷ content of the plasma. The emissivity ratio for the DR/3.12 keV has its maximum value of 0.04 at kT = 0.7 keV, but the emissivity of both lines is weak here, so any hotter component will dominate and lead to a lower ratio being observed.

To avoid cool gas in the Perseus core contaminating the flux of the nearby Ar and K lines, we also tried excising the central region with 1' radius of the cluster and performed the fit on the core-excised co-added MOS spectrum. We found that adding an extra Gaussian line at 3.57 keV has improved the fit by Δχ² of 12.8 for an additional degree of freedom with a best-fit flux of 2.1 $^{+0.7}_{-0.6}$ ($^{+1.2}_{-1.1}$) × 10⁻⁵ photons cm⁻² s⁻¹ (see Figure 8). Excising the innermost 1' reduced the flux of the detected line by a factor of two, indicating that most of the flux of this emission originates from the cool core. The mixing angle that corresponds to the line flux from the core-excised Perseus spectrum is consistent within 1σ–2σ with those for the bright clusters (Centaurus+Coma+Ophiuchus) and the full sample, respectively (Table 5).

Zoom In Zoom Out Reset image size

We also note that some scatter of the dark matter decay signal between individual clusters is naturally expected. For example, one can imagine a filament of dark matter along the line of sight in the direction of Perseus, which may boost the flux of the detected line and cause tension between the Perseus Cluster and the full sample. However, such a filament would have to be rather extreme.

With the knowledge that the 3.62 keV line can be anomalously high (at least in Perseus), we should now try to re-fit the stacked MOS spectrum of the full sample to see if the line in the full sample is affected by the 3.62 keV excess from Perseus, which is part of the full sample. We set the flux of the 3.62 keV line to the Perseus contribution of the Ar xvii DR line to the full-sample spectrum (2.3 × 10⁻⁶ photons cm⁻² s⁻¹), assuming that all the new line flux in Perseus originates from the abnormally bright DR line. We note that this flux was already a factor of 30 above the predicted upper limits by AtomDB. Adding an extra Gaussian component, representing the new line, to a model with the anomalous 3.62 keV line still improves the fit by Δχ² of 6.52 for 2 dof. The best-fit energy and flux were 3.55 ± 0.03 (0.05) and 2.2$^{+1.6}_{-0.9}$ ($^{+2.2}_{-1.5}$) × 10⁻⁵ photons cm⁻² s⁻¹, respectively. The new line is still required with 2.5σ in the full sample; however, the energy of this line gets lower and its confidence interval wider. The line energy comes into agreement with the energy detected in the PN full sample (see Figure 9, left panel). If we completely free the normalization of the 3.62 keV line in the full-sample MOS spectrum, it becomes lower than the Perseus contribution that we considered above.

Zoom In Zoom Out Reset image size

Each spectrum was fitted using a standard multi-temperature apec model as described in Section 2.3 to determine the best-fit X-ray redshift of each observation, shown in Table 6. Each event file was then blueshifted to the cluster's source frame using these best-fit redshifts. The source and background spectra in the source's frame were obtained by generating spectra using the scaled event energy values in the event files. The ARFs and RMFs were remapped based on the estimates of the best-fit redshifts. The RMFs and ARFs were weighted by only the exposure time of each observation. The scaled source and background spectra were co-added using the FTOOL mathpha, whereas ARFs and RMFs were merged using the FTOOLS addarf and addrmf tools, respectively.

Following the same method as described in Section 3.1, the continuum emission was fit using the line-free apec model with additional Gaussian models to represent the strong emission lines. The best-fit temperature from the 2.0–6.0 keV band and normalizations of the line-free apec model, fluxes, and equivalent widths of S xvi, Ca xix, and Ca xx lines at 2.63 keV, 3.90 keV, and 4.11 keV are given in Table 7. We have searched especially the 3.0–4.0 keV interval where the 3.57 keV line emission was detected in the stacked XMM-Newton observations. The measured fluxes of S xvi, Ca xix, and Ca xx lines from the ACIS-I and ACIS-S spectra with the AtomDB fluxes yielded the maximum predicted fluxes of K xviii lines at 3.47 keV and 3.51 keV, Ar xvii line at 3.68 keV, and K xviii line at 3.71 keV as described in detail in Section 3.1. The triplet emission line at Ar xvii 3.12 keV was used to determine the maximum allowed flux of the Ar xvii DR line at 3.62 keV at any plasma temperature as described above. The predicted fluxes of these lines are given in Table 8. Using 0.1 and 3 times the upper bound of these estimates as lower and upper limits for K xviii and Ar xvii and 10⁻³–10⁻² times the flux of the Ar xvii triplet for the lower and upper bounds for the Ar xvii DR line, we determined the best-fit flux of the weak residual around 3.57 keV.

An additional Gaussian model improves the fit by Δχ² of 11.8 for an additional 2 dof. The line was unresolved and consistent with broadening by the instrument response in the Perseus Cluster spectra. The Perseus ACIS-S spectra yield a best-fit energy of 3.56 ± 0.02 (0.03) keV for an additional Gaussian model, given in Table 5. The flux of the detected signal is 1.02$^{+0.4}_{-0.4}$ ($^{+0.5}_{-0.5}$) × 10⁻⁵ $\rm {photons}\ \rm {cm}^{-2}\ \rm {s}^{-1}$. This detection corresponds to a false detection probability of 0.5% in the co-added ACIS-S spectrum. The right panel of Figure 10 shows the signal in the Chandra ACIS-S observations of the Perseus Cluster before and after the Gaussian model is added to the fit.

Zoom In Zoom Out Reset image size

To further demonstrate that the detected flux is independent of the spectral modeling, we fit the ACIS-S spectrum of the Perseus Cluster with a two-temperature vapec model with abundances of trace elements set to that of Fe. We obtained an acceptable fit in the 3–6 keV energy band with χ² of 182.1 for 147 dof. An additional Gaussian model at 3.56 keV (rest energy) improved the fit by Δχ² of 16 for an extra degree of freedom. The best-fit flux of the line is 1.09 ± 0.26 (0.42) × 10⁻⁵ photons cm² s⁻¹, which is consistent with the flux measured in the line-free apec model fit with additional Gaussian models. This test shows that the detection is robust and independent of the method used in the spectrum modeling. The Perseus co-added spectrum fit with a two-temperature vapec model is shown in Figure 11.

Zoom In Zoom Out Reset image size

We then performed the same search in the co-added ACIS-I spectrum of the Perseus Cluster. Fitting the 2.5–6 keV band of the ACIS-I spectrum with a line-free apec model with additional Gaussian lines as described above produced a good fit overall with a total χ² of 158.7 for 152 dof. Adding a Gaussian line at 3.56 keV, the energy where the line was detected in the co-added ACIS-S spectrum of the Perseus Cluster, improved the fit by Δχ² of 6.2 for an additional degree of freedom. The flux of the detected signal was 1.9 $^{+0.8}_{-0.8}$ ($^{+1.2}_{-1.6}$) × 10⁻⁵ $\rm {photons}\ \rm {cm}^{-2}\ \rm {s}^{-1}$ in the co-added ACIS-I spectrum. The left panel of Figure 10 shows the ACIS-I spectrum of the Perseus Cluster before and after an additional Gaussian model is added to the total model to demonstrate the detection of the line.

The mixing angle sin ²(2θ) estimate from the co-added Chandra ACIS-S observations of the Perseus Cluster is 4.0 $^{+1.5}_{-1.4}$ ($^{+1.8}_{-1.8}$) × 10⁻¹⁰, which is consistent with the angle obtained from the co-added ACIS-I and XMM-Newton MOS observations of the Perseus Cluster at the 1σ level. Since the ACIS-S chip covers the central 4' region of the Perseus core, higher flux measured from ACIS-S observations also indicates that this emission is concentrated in the core, confirming the results from the XMM-Newton observations of the Perseus core.

We have performed the same fitting strategy described above to the co-added spectra of the Virgo Cluster, e.g., line-free apec model with additional Gaussian lines. We used the lower and upper limits to the K and Ar line in the 3.4–3.7 keV band based on the upper limits estimated from AtomDB (given in Table 8). The overall fit was acceptable with a total χ² of 82.5 for 62 dof. Unlike the Perseus Cluster, the co-added Virgo Cluster did not show any residuals around 3.57 keV in the fit with the line-free apec model. Adding a Gaussian line did not significantly improve the fit. We were able to place an upper limit of 9.1 × 10⁻⁶$\rm {photons}\ \rm {cm}^{-2}\ \rm {s}^{-1}$ at the 90% confidence level. This limit corresponds to an upper limit on the mixing angle of sin ²(2θ) < 1.1 × 10⁻¹⁰.

We also fit the 2.5–4.0 keV band of the Virgo spectrum using a two-temperature standard vapec model. The fit has a total χ² obtained from the vapec model of 91.7 for 82 dof. We overall obtained a better fit with the standard vapec model than the fit with the line-free apec model. The best-fit model also did not require the addition of a line at 3.56 keV. The 90% upper limit to the flux of this line is <6.2 × 10⁻⁶ photons cm⁻² s⁻¹. The differences in the modeling approaches used in the ACIS-I spectrum fits of the Virgo Cluster (line-free apec with Gaussians and vapec) are demonstrated in Figure 12. The factor of two difference in the upper limits on the flux measurements indicates that the systematical uncertainties in the flux measurements can be as large as a factor of two depending on the modeling method used in this analysis.

Zoom In Zoom Out Reset image size

One possible interpretation is that the detected line is an unknown plasma emission line. The flux of the line corresponds to a maximum emissivity of 3.3 × 10⁻¹⁸ photons cm³ s⁻¹, derived using the emission measure appropriate for the lowest temperature (4.36 keV) component as described in Section 3.1. For comparison, this is similar to the maximum emissivity of the Ca xx Lyα line at 4.1 keV. Given that the Ca xx line was previously observed in individual galaxy cluster spectra, including the Perseus Cluster (e.g., Tamura et al. 2009), a line as strong at ∼3.57 keV would have been observed had it been expected. However, there is no likely candidate for an atomic transition near 3.57 keV. The emission lines of strong hydrogen- and helium-like ions are well known, and none fall in this band. The only candidate emission line from such an ion would be the He-like Cl xvi n = 5 → 1 transition at 3.52 keV, but if this is the cause, it would imply the presence of even stronger lines from the n = 3 → 1 and n = 4 → 1 transitions at 3.27 and 3.44 keV, respectively, but these are not seen. Emission lines from L-shell ions form a far more complex pattern. However, the binding energy of Li-like Zn (Z = 30) is only 2.782 keV, so the transition lines of all lighter elements or less ionized species must be at lower energies than this. If this line is a K-shell fluorescence transition, it must be from an element whose neutral and Li-like K-shell fluorescent line energies bound 3.57 keV. The only such atoms are argon and potassium, but in this case the relevant Ar K-shell fluorescence transition is simply another name for the Ar xvii DR line discussed in detail above. The neutral potassium Kα fluorescence line is at 3.313 keV, while neutral Kβ is at 3.59 keV, so there must be transitions at the relevant energy. In this case, the best matches are the Kα transitions of K xvi through K xiv ions, which occur at ∼3.57 keV (Palmeri et al. 2012). However, since at any temperature above 1 keV potassium will have at most two bound electrons, any such line would have to be originating from an unknown source of photoionized potassium in clusters. Thus, this scenario is very unlikely, since the compact sources (e.g., active galactic nuclei) are not strong enough to photoionize the low-density ICM.

Although a complete analysis was not shown, adding an Ar xvii DR line at 3.62 keV with unconstrained flux into all of our spectra would significantly impact both the fit results and detection level of a line at 3.57 keV. We have constrained this line to be at most 1% of the strength of the unresolved Ar xvii triplet at 3.12 keV, but must consider the physical situation required to maximize the 3.62 keV DR line. In thermal equilibrium, the maximum strength of this line is 4% of the Ar xvii triplet, albeit at a temperature where the expected emission is negligible. One might also consider an extreme nonequilibrium situation with cold electrons that are unable to collisionally excite any Ar xvii lines, but DR is still possible. Examining the satellite line data in AtomDB, taken from Vainshtein & Safronova (1980), shows that even in this case the maximum ratio is only 7%, as there are DR satellite lines at the energies of the Ar xvii triplet as well and these lines would also be excited in such a case. While not physically impossible if there was a significant and unexpected error in the atomic physics calculations, we have no reason to believe that this has occurred.

One other possibility is a radiative recombination continuum (RRC) edge feature. The S xvi recombination edge lies at 3.494 keV, and if it was bright enough, it might fill in some of the flux in this region. However, we note that the use of the no-line model has already included the RRC feature under equilibrium conditions. Producing a stronger RRC would require a sharp drop in the electron temperature, while retaining a large, fully stripped S^{16 +} population. Also, if the temperature drops below ≈0.1 keV, the RRC feature becomes very narrow and will be an order of magnitude less powerful at 3.57 keV compared with right at the edge: this shape is not consistent with our observations. Similarly, at hotter temperatures, the RRC becomes almost constant with energy once above the edge. If the RRC was really there, we would expect to see a residual at about 3.6 keV, which we do not. Finally, we note that the edge is 50–80 eV from our proposed line, which makes an unlikely source of the line.

We also note that our assumptions regarding relative line strengths have assumed that the ICM is in thermal equilibrium or close to it. Charge exchange between highly ionized ions and neutral hydrogen or helium could also create X-ray emission lines with different ratios (Smith et al. 2012). This could affect our assumption of equilibrium line ratios, although we have included a substantial range around the equilibrium values. It is important to note that these CX lines are not new,but rather the same lines occurring in different ratios. Due to its large cross section relative to electron excitation rates, astrophysical CX can occur only in a thin sheet where ions and neutrals interact directly, limiting its total emission relative to the large ICM volume. In certain cases, such as the core of the Perseus Cluster where many neutral filaments are known, it is possible that CX could be large enough to create a small fraction of the total X-ray emission, although it would not create or enhance a line at 3.57 keV or the DR line at 3.62 keV. CX could not dominate the overall emission, however, as it would also create Fe xvii and other lines that are not detected.

An interesting interpretation of the line is the decay signature of the sterile neutrino, a long-sought dark matter particle candidate (e.g., Boyarsky et al. 2009; see our Section 1). The mass of the sterile neutrino would be double the decay photon energy, m_s = 7.1 keV. The line flux detected in our full sample corresponds to a mixing angle for the decay sin ²(2θ) ∼ 7 × 10⁻¹¹. This value is below the upper limits placed by the previous searches, shown in Figure 13. Our detections from the stacked XMM-Newton MOS observations of galaxy clusters are shown with a star in red in that figure. Figure 14 shows the detections and upper limits we obtained from our various subsamples we used in this work (based on the included cluster masses and distances), as well as a comparison with the previous upper limit placed using the Bullet Cluster by Boyarsky et al. (2008) at 3.57 keV, which is the most relevant earlier constraint for us. Since the mixing angle is a universal quantity, all the subsample measurements must agree.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The line in the subsample of 69 fainter clusters (full sample sans Perseus, Coma, Ophiuchus, and Centaurus) corresponds to a mixing angle that is consistent with the full sample; the same is seen (though with a mild 1.5σ tension) for the subsample of bright nearby clusters Coma+Centaurus+Ophiuchus. However, the brightness of the new line in the XMM-Newton spectrum of Perseus corresponds to a significantly higher mixing angle than that for the full sample (by factor of eight in terms of the line flux for the MOS spectrum), which poses a problem in need of further investigation. In principle, an enhanced flux of the detected line in the Perseus spectra may be due to a dark matter filament along the line of sight, though it would have to be rather extreme, so it is wise to look for more plausible explanations.

We tried to excise the central 1' region of the Perseus Cluster, to see if the flux originates in the cool core of the cluster. Indeed, this decreased the flux in the line in half and removed most of the tension with the other measurements. However, this suggests that either some of the line flux is astrophysical in origin (at least in Perseus) or the cool gas in the core of the cluster affects our ability to measure the continuum and the fluxes of the nearby K xviii and Ar xvii lines, in the end resulting in an overestimate of the flux of our detected line. It appears that in Perseus, there is an anomalously strong line at the position of the Ar xvii DR line at 3.62 keV.

With this knowledge, we have tried to add this anomalous 3.62 keV line in the model for the full sample, where we have the most statistically significant line detection. The additional line is still required, albeit at a lower significance and a slightly lower energy of 3.55 ± 0.03 keV. Note that the sample of bright clusters is dominated by the emission from the cool cores of the Ophiuchus and Centaurus Clusters; if this Ar 3.62 keV line anomaly is typical of cool cores, they may also be affected. However, freeing the flux of the 3.62 keV line in the MOS full-sample fit did not require additional contribution from clusters other than Perseus, though the constraints are obviously weak.

The radial distribution of the flux of this line should be investigated further in the nearby bright clusters, including those with and without cool cores.

We note that even if the sterile neutrino interpretation of the emission line is correct, this detection would not necessarily imply that all dark matter is composed of these particles. Assuming a standard cosmological history below a temperature of a few hundred MeV, sterile neutrinos would be produced by oscillations with active neutrinos at an abundance determined by the mass and mixing angle (e.g., Dodelson & Widrow 1994; Kusenko 2009). Accounting for the increase in mixing angle that would be inferred for a dark matter fraction in sterile neutrinos less than unity, we find that this fraction is ∼13%−19% based on the methods in Abazajian (2006) and Asaka et al. (2007)—and cannot exceed 26% based on the absolute lower bound distorted wave production estimate in Asaka et al. (2007).

This implies that either (1) sterile neutrinos are a subdominant component of dark matter, (2) sterile neutrinos are predominantly produced by some other mechanism, or (3) the emission line originates from some other radiatively decaying light dark matter candidate such as moduli dark matter (Kusenko et al. 2013). The Shi−Fuller mechanism is one of the possible production mechanisms for the sterile neutrino dark matter interpretation of this detection. The implications of the detection for structure formation in cosmological small scales are discussed in detail in Abazajian (2014).

They may also be produced by means that do not involve oscillations, such as inflaton or Higgs decay (Kusenko 2006, 2009; Shaposhnikov & Tkachev 2006; Petraki & Kusenko 2008), although there may still be sufficient mixing to provide an observable radiative decay signal. This detection is consistent with 100% of dark matter composed of sterile neutrinos produced by these mechanisms, as well as by the split seesaw mechanism (Kusenko et al. 2010). Even in this case, some sterile neutrinos would be produced by nonresonant oscillations. However, based again on the calculations in Abazajian (2006) and Asaka et al. (2007), only 1%−3% of the sterile neutrino abundance (with an upper limit of 7%) would be accounted for in this way for a sterile neutrino with mass of 7.1 keV and a mixing angle corresponding to sin ²(2θ) ∼ 7 × 10⁻¹¹.

Our result must be verified using a variety of X-ray instruments, X-ray-emitting dark-matter-dominated objects, methods of data reduction, background subtraction, and statistical techniques to investigate the interpretation of this line. The future high-resolution Astro-H observations will be able to measure the broadening of the line, which will allow us to measure its velocity dispersion. To detect a dark matter decay line, which is much weaker than the plasma lines, will require a significantly long exposure. We performed 1 Ms Astro-H soft X-ray spectrometer (SXS) simulations of the Perseus Cluster assuming that the width (15 eV) of the dark matter decay line is determined by the virial velocities of dark matter particles of 1300 km s⁻¹. Figure 15 shows that the broader dark matter line will be easily distinguished from the plasma emission lines, which are only broadened by the turbulence in the X-ray-emitting gas.

Zoom In Zoom Out Reset image size