Dark Matter Strikes Back at the Galactic Center

TL;DR

This paper investigates how unmodeled astrophysical sources can bias the interpretation of the Galactic Center Excess, revealing that dark matter signals might be underestimated due to modeling limitations.

Contribution

It demonstrates that unmodeled sources can cause misattribution of dark matter signals to point sources in non-Poissonian template fitting analyses.

Findings

Unmodeled sources in simulations can mislead the identification of dark matter signals.

Real data shows effects consistent with source mis-modeling, affecting dark matter attribution.

Dark matter may still be the dominant contributor to the Galactic Center Excess.

Abstract

Statistical evidence has previously suggested that the Galactic Center GeV Excess (GCE) originates largely from point sources, and not from annihilating dark matter. We examine the impact of unmodeled source populations on identifying the true origin of the GCE using non-Poissonian template fitting (NPTF) methods. In a proof-of-principle example with simulated data, we discover that unmodeled sources in the Fermi Bubbles can lead to a dark matter signal being misattributed to point sources by the NPTF. We discover striking behavior consistent with a mismodeling effect in the real Fermi data, finding that large artificial injected dark matter signals are completely misattributed to point sources. Consequently, we conclude that dark matter may provide a dominant contribution to the GCE after all.

Tables

Table 1. Table 1: Bayes factors and 95% containment intervals on the DM flux posterior, for various analyses on simulated data. Note that in addition to the templates listed, in all cases we also simulate smooth Bubbles, Isotropic, and Galactic Diffuse templates. In the simulation, the DM flux is not correctly recovered in the case where the fit is performed with the standard pipeline (shown in purple). The data analysis that is favored has the relative Bayes factor appear next to the analysis row, where the two analyses that are being compared are not grayed out.

Simulated Data, 3FGL Masked
Simulation	Simulated	Analysis Templates	DM Flux	Bayes Factor
	DM Flux		(95%)
Bubbles PS		Same as simulated	$[\mathrm{\hspace{0.17em}1.3},2.5]\%$	$7\times {10}^{20}$		$1\times {10}^{26}$
Disk PS	$2\%$	Same but Bubbles PS $\to$ NFW PS	$\left[\text{0.0, 0.4}\right]\%$		$2\times {10}^5$
DM		Same but no Bubbles PS	$[\mathrm{\hspace{0.17em}0.8},2.2]\%$
Bubbles PS		Same as simulated	$[\mathrm{\hspace{0.17em}4.6},5.8]\%$	$9\times {10}^{16}$		$8\times {10}^{23}$
Disk PS	$5\%$	Same but Bubbles PS $\to$ NFW PS	$\left[\text{0.0, 1.8}\right]\%$		$9\times {10}^6$
DM		Same but no Bubbles PS	$[\mathrm{\hspace{0.17em}4.3},5.4]\%$
Bubbles PS		Same as simulated	$[\mathrm{\hspace{0.17em}14.6},15.8]\%$	$5\times {10}^{16}$		$5\times {10}^{21}$
Disk PS	$15\%$	Same but Bubbles PS $\to$ NFW PS	$\left[\text{11.4, 13.7}\right]\%$		$1\times {10}^5$
DM		Same but no Bubbles PS	$[\mathrm{\hspace{0.17em}14.3},15.6]\%$

Table 2. Table 2: Bayes factors for analyses of the real Fermi data injected with a DM signal. The DM flux is not correctly recovered when the standard PS templates are considered (shown as purple), suggesting a large downward bias to the DM flux that could hide a real signal, similar to the proof-of-principle case detailed in Tab.  1 . The data analysis that is favored has the relative Bayes factor appear next to the analysis row, where the two analyses that are being compared are not grayed out.

Real Data, 3FGL Masked
Injected	Analysis Templates	DM Flux	Bayes Factor
DM Flux		(95%)
None	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}1.1},2.4]\%$
	Diffuse + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.3}\right]\%$	$1\times {10}^9$
	Diffuse + Iso + Bub + DM
$1.8\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}2.8},4.2]\%$
	Diffuse + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.3}\right]\%$	$2\times {10}^{12}$	$4\times {10}^2$
	Diffuse + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	Fixed at injection			$7\times {10}^9$
	Diffuse + Iso + Bub + Fixed DM	value ($1.8\%$)
$6.7\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}7.7},9.0]\%$
	Diffuse + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 1.3}\right]\%$	$7\times {10}^{11}$
	Diffuse + Iso + Bub + DM
$15.2\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}15.7},16.9]\%$
	Diffuse + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{5.0, 11.5}\right]\%$	$8\times {10}^9$
	Diffuse + Iso + Bub + DM

Table 3. Table S1: Parameters and associated prior ranges for the Inner Galaxy analyses, when using the Fermi p6v11 diffuse model, as in the main text.

Prior Ranges
Parameter	Inner Galaxy (masked)	Inner Galaxy (unmasked)
${\log }_{10}{A}_{\text{iso}}$	$[-3,1]$	$[-3,1]$
${\log }_{10}{A}_{\text{dif}}$	$[0,2]$	$[0,2]$
${\log }_{10}{A}_{\text{bub}}$	$[-3,1]$	$[-3,1]$
${\log }_{10}{A}_{\text{NFW}}$	$[-3,1]$	$[-3,1]$
${\log }_{10}{A}_{\text{PS}}^{\text{NFW}}$	$[-6,1]$	$[-6,1]$
$S_b^{\text{NFW}}$	$[0.05,60]$	$[0.05,60]$
$n_1^{\text{NFW}}$	$[2.05,60]$	$[5,45]$
$n_2^{\text{NFW}}$	$[-3,1.95]$	$[-3,1.95]$
${\log }_{10}{A}_{\text{PS}}^{\text{disk}}$	$[-6,1]$	$[-6,1]$
$S_b^{\text{disk}}$	$[0.05,60]$	$[0.05,60]$
$n_1^{\text{disk}}$	$[2.05,60]$	$[2.05,5]$
$n_2^{\text{disk}}$	$[-3,1.95]$	$[-3,1.95]$
${\log }_{10}{A}_{\text{PS}}^{\text{iso}}$	$[-6,1]$	$[-6,1]$
$S_{b_1}^{\text{iso}}$	$[5,40]$	$[1,40]$
$S_{b_2}^{\text{iso}}$	$[0.05,30]$	$[0.05,30]$
$n_1^{\text{iso}}$	$[2.05,5]$	$[2.05,5]$
$n_2^{\text{iso}}$	$[1.5,4.5]$	$[0.5,4.5]$
$n_3^{\text{iso}}$	$[-1.95,1.95]$	$[-1.95,1.95]$

Table 4. Table S2: Parameter values used to generate the simulated data from Poissonian and non-Poissonian templates for the proof-of-principle scenario, for cases where the 3FGL PSs are masked or unmasked. Note that for simplicity, in the proof-of-principle scenario, the Isotropic PS are not simulated when the 3FGL are masked, however they are included in the unmasked case to ensure the bright 3FGL are correctly captured.

Simulation Parameters
Parameter	Inner Galaxy (masked)	Inner Galaxy (unmasked)
${\log }_{10}{A}_{\text{iso}}$	$-1.94$	$-2.20$
${\log }_{10}{A}_{\text{dif}}$	$1.11$	$1.11$
${\log }_{10}{A}_{\text{bub}}$	$-0.08$	$-0.06$
${\log }_{10}{A}_{\text{NFW}}$	$-0.56$	$-0.42$
${\log }_{10}{A}_{\text{PS}}^{\text{NFW}}$, ${\log }_{10}{A}_{\text{PS}}^{\text{bub}}$	$-2.08$	$-2.37$
$S_b^{\text{NFW}}$, $S_b^{\text{bub}}$	13.47	$20.53$
$n_1^{\text{NFW}}$, $n_1^{\text{bub}}$	35.13	$32.98$
$n_2^{\text{NFW}}$, $n_2^{\text{bub}}$	$-1.68$	$-1.41$
${\log }_{10}{A}_{\text{PS}}^{\text{disk}}$	$-4.15$	$-3.71$
$S_b^{\text{disk}}$	37.69	$42.70$
$n_1^{\text{disk}}$	31.59	$2.40$
$n_2^{\text{disk}}$	$-0.88$	$-1.50$
${\log }_{10}{A}_{\text{PS}}^{\text{iso}}$	$-4.78$	$-4.18$
$S_{b,1}^{\text{iso}}$	$22.59$	$23.85$
$S_{b,2}^{\text{iso}}$	$3.22$	$5.22$
$n_1^{\text{iso}}$	$3.89$	$3.41$
$n_2^{\text{iso}}$	$3.84$	$3.48$
$n_3^{\text{iso}}$	$-0.27$	$-0.18$

Table 5. Table S3: Same as the simulated data proof-of-principle scenario in Tab.  1 , but an unmasked analysis. Purple signals the standard pipeline analysis.

Simulated Data, 3FGL Unmasked
Simulation	Injected	Analysis Templates	DM Flux	Bayes Factor
	DM Flux		(95%)
Bubbles PS		Same as simulated	$[\mathrm{\hspace{0.17em}2.3},3.1]\%$	$1\times {10}^{34}$		$4\times {10}^{36}$
Disk PS	$2.4\%$	Same but Bubbles PS $\to$ NFW PS	$\left[\text{0.0, 1.7}\right]\%$		$3\times {10}^2$
DM		Same but no Bubbles PS	$[\mathrm{\hspace{0.17em}2.4},3.4]\%$

Table 6. Table S4: Bayes factors for analyses of the real Fermi data injected with a DM signal. The DM flux is not correctly recovered when the standard PS templates are considered (shown as purple), similar to the proof-of-principle case detailed in Tabs.  1  and  S3 . The data analysis that is favored has the relative Bayes factor appear next to the analysis row, where the two analyses that are being compared are not grayed out.

Real Data, 3FGL Unmasked
Injected	Analysis Templates	DM Flux	Bayes Factor
DM Flux		(95%)
None	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}1.9},2.9]\%$
	Diffuse + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.2}\right]\%$	$4\times {10}^{18}$
	Diffuse + Iso + Bub + DM
$2.4\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}3.7},4.7]\%$
	Diffuse + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.3}\right]\%$	$6\times {10}^{19}$	$9\times {10}^2$
	Diffuse + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	Fixed at injection			$6\times {10}^{16}$
	Diffuse + Iso + Bub + Fixed DM	value ($2.4\%$)

Table 7. Table S5: Bayes factors for analyses of simulated data injected with a DM signal. This is the simulated data parallel of Tab.  2 , assuming the best-fit contains no DM, and injecting DM on top of the best-fit. In this case, the DM flux shows distinctly different behavior to the real data case when the standard PS templates are considered (shown in purple); while it is still not always correctly reconstructed, the discrepancy is much smaller than in the real data. The data analysis that is favored has the relative Bayes factor appear next to the analysis row, where the two analyses that are being compared are not grayed out.

Simulated Data, 3FGL Masked
Simulation	Injected	Analysis Templates	DM Flux	Bayes Factor
	DM Flux		(95%)
	None	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}0.4},1.9]\%$
		Diffuse + Iso + Bub + DM
		Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.3}\right]\%$	$8\times {10}^6$
		Diffuse + Iso + Bub + DM
	$1.8\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}2.4},3.9]\%$
Disk PS		Diffuse + Iso + Bub + DM
Iso PS		Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 2.1}\right]\%$	$1\times {10}^6$	$1$
NFW PS		Diffuse + Iso + Bub + DM
Diffuse		Disk PS + Iso PS + NFW PS	Fixed at injection			$1\times {10}^6$
Iso		Diffuse + Iso + Bub + Fixed DM	value ($1.8\%$)
Bub	$6.7\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}7.8},9.2]\%$
		Diffuse + Iso + Bub + DM
		Disk PS + Iso PS + NFW PS	$\left[\text{5.0,7.6}\right]\%$	$8\times {10}^5$
		Diffuse + Iso + Bub + DM
	$15.2\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}12.8},15.4]\%$
		Diffuse + Iso + Bub + DM
		Disk PS + Iso PS + NFW PS	$\left[\text{10.0, 16.5}\right]\%$	$1\times {10}^5$
		Diffuse + Iso + Bub + DM

Table 8. Table S6: Parameters and associated prior ranges for the Inner Galaxy analyses, when using diffuse models A and F, where the inverse Compton and π 0 + limit-from superscript 𝜋 0 \pi^{0}+ bremsstrahlung contributions are considered independently with their own templates.

Prior Ranges, Inner Galaxy, Diffuse Models A and F
Parameter	Model A (masked)	Model A (unmasked)	Model F (masked)	Model F (unmasked)
${\log }_{10}{A}_{\text{iso}}$	$[-3,1]$	$[-3,1]$	$[-3,1]$	$[-3,1]$
${\log }_{10}{A}_{\text{dif}}^{\text{ICS}}$	$[-2,1]$	$[-2,1]$	$[-2,1]$	$[-2,1]$
${\log }_{10}{A}_{\text{dif}}^{{\pi }^0\text{brem}}$	$[-1,2]$	$[-1,2]$	$[-1,2]$	$[-1,2]$
${\log }_{10}{A}_{\text{bub}}$	$[-3,1]$	$[-3,1]$	$[-3,1]$	$[-3,1]$
${\log }_{10}{A}_{\text{NFW}}$	$[-3,1]$	$[-3,1]$	$[-3,1]$	$[-3,1]$
${\log }_{10}{A}_{\text{PS}}^{\text{NFW}}$	$[-6,1]$	$[-6,1]$	$[-6,1]$	$[-6,1]$
$S_b^{\text{NFW}}$	$[0.05,60]$	$[0.05,30]$	$[5,60]$	$[0.05,60]$
$n_1^{\text{NFW}}$	$[50,95]$	$[5,45]$	$[2.05,60]$	$[5,60]$
$n_2^{\text{NFW}}$	$[-3,1.95]$	$[-3,1.95]$	$[-3,1.95]$	$[-3,1.95]$
${\log }_{10}{A}_{\text{PS}}^{\text{disk}}$	$[-6,1]$	$[-6,1]$	$[-6,1]$	$[-6,1]$
$S_b^{\text{disk}}$	$[0.05,60]$	$[0.05,50]$	$[0.05,60]$	$[0.05,60]$
$n_1^{\text{disk}}$	$[25,80]$	$[2.05,5]$	$[25,80]$	$[2.05,65]$
$n_2^{\text{disk}}$	$[-3,1.95]$	$[-3,1.95]$	$[-3,1.95]$	$[-3,1.95]$
${\log }_{10}{A}_{\text{PS}}^{\text{iso}}$	$[-6,1]$	$[-6,1]$	$[-6,1]$	$[-6,1]$
$S_{b_1}^{\text{iso}}$	$[10,50]$	$[1,40]$	$[1,40]$	$[1,40]$
$S_{b_2}^{\text{iso}}$	$[0.05,30]$	$[0.05,30]$	$[0.05,30]$	$[0.05,30]$
$n_1^{\text{iso}}$	$[2.05,5]$	$[2.05,5]$	$[2.05,5]$	$[2.05,5]$
$n_2^{\text{iso}}$	$[1.5,4.5]$	$[0.5,4.5]$	$[0.5,4.5]$	$[0.5,4.5]$
$n_3^{\text{iso}}$	$[-1.95,1.95]$	$[-1.95,1.95]$	$[-1.95,1.95]$	$[-1.95,1.95]$

Table 9. Table S7: Bayes factors for analyses of the real Fermi data injected with a DM signal, in the Inner Galaxy ROI (unmasked), using the diffuse model A. The standard-pipeline analysis (modulo the replacement of the diffuse model) is shown in purple. The data analysis that is favored has the relative Bayes factor appear next to the analysis row, where the two analyses that are being compared are not grayed out.

Diffuse Model A, Real Data, 3FGL Unmasked
Injected	Analysis Templates	DM Flux	Bayes Factor
DM Flux		(95%)
None	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}4.0},5.0]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.7}\right]\%$	$3\times {10}^9$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub+ DM
$2\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}5.7},6.8]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.8}\right]\%$	$3\times {10}^9$	$10$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	Fixed at injection			$5\times {10}^8$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + Fixed DM	value ($2\%$)

Table 10. Table S8: Bayes factors for analyses of the real Fermi data injected with a DM signal, in the Inner Galaxy ROI (unmasked), using the diffuse model F. The standard-pipeline analysis (modulo the replacement of the diffuse model) is shown in purple. The data analysis that is favored has the relative Bayes factor appear next to the analysis row, where the two analyses that are being compared are not grayed out.

Diffuse Model F, Real Data, 3FGL Unmasked
Injected	Analysis Templates	DM Flux	Bayes Factor
DM Flux		(95%)
None	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}2.6},3.7]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.6}\right]\%$	$6\times {10}^3$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
$2\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}4.2},5.8]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.9}\right]\%$	$1\times {10}^5$	$10$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	Fixed at injection			$1\times {10}^4$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + Fixed DM	value ($2\%$)

Table 11. Table S9: Bayes factors for analyses of the real Fermi data injected with a DM signal, in the Inner Galaxy ROI (masked), using the diffuse model A. The standard-pipeline analysis (modulo the replacement of the diffuse model) is shown in purple. The data analysis that is favored has the relative Bayes factor appear next to the analysis row, where the two analyses that are being compared are not grayed out.

Diffuse Model A, Real Data, 3FGL Masked
Injected	Analysis Templates	DM Flux	Bayes Factor
DM Flux		(95%)
None	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}4.0},5.0]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.1}\right]\%$	$1\times {10}^3$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
$2\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}5.5},6.7]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 1.0}\right]\%$	$2\times {10}^3$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM

Table 12. Table S10: Bayes factors for analyses of the simulated data (based on the best fit with NFW PS explaining the GCE) injected with a DM signal, in the Inner Galaxy ROI (masked), using the diffuse model A. The standard-pipeline analysis (modulo the replacement of the diffuse model) is shown in purple. The data analysis that is favored has the relative Bayes factor appear next to the analysis row, where the two analyses that are being compared are not grayed out.

Diffuse Model A, Simulated Data, 3FGL Masked
Injected	Analysis Templates	DM Flux	Bayes Factor
DM Flux		(95%)
None	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}4.2},5.6]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.4}\right]\%$	$5\times {10}^4$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
$2\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}6.0},7.5]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 0.5}\right]\%$	$2\times {10}^4$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM

Table 13. Table S11: Bayes factors for analyses of the real Fermi data injected with a DM signal, in the Inner Galaxy ROI (masked), using the diffuse model F. The standard-pipeline analysis (modulo the replacement of the diffuse model) is shown in purple.The data analysis that is favored has the relative Bayes factor appear next to the analysis row, where the two analyses that are being compared are not grayed out.

Diffuse Model F, Real Data, 3FGL Masked
Injected	Analysis Templates	DM Flux	Bayes Factor
DM Flux		(95%)
None	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}1.2},3.0]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 1.0}\right]\%$	$2$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
$2\%$	Disk PS + Iso PS	$[\mathrm{\hspace{0.17em}3.5},4.7]\%$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM
	Disk PS + Iso PS + NFW PS	$\left[\text{0.0, 1.2}\right]\%$	$1$
	Diff IC + Diff ${\pi }^0$brem + Iso + Bub + DM

Equations

$${\mu }_{p,\ell }(\mathbf{\theta })=A_{\ell }(\mathbf{\theta })T_{p,\ell },$$

$${\mu }_p^{n_p}(\mathbf{\theta })=\sum _{\ell }{\mu }_{p,\ell }(\mathbf{\theta }),$$

$$P_{n_p}^{(p)}(\mathbf{\theta })=\frac{{\mu }_p^{n_p}(\mathbf{\theta })}{n_p!}{e}^{-{\mu }_p(\mathbf{\theta })}.$$

$$P(d|\mathbf{\theta },\mathcal{M})=\prod _p{P}_{n_p}^{(p)}(\mathbf{\theta }).$$

$$\frac{dN_{p}}{dF}(F;{\bm{\theta}})=A({\bm{\theta}})T^{({\rm PS})}_{p}\left\{\begin{array}[]{lc}\left(\frac{F}{F_{b,1}}\right)^{-n_{1}},&F\geq F_{b,1}\\ \left(\frac{F}{F_{b,1}}\right)^{-n_{2}},&F_{b,1}>F\geq F_{b,2}\\ \left(\frac{F_{b,2}}{F_{b,1}}\right)^{-n_{2}}\left(\frac{F}{F_{b,2}}\right)^{-n_{3}},&F_{b,2}>F\geq F_{b,3}\\ \left(\frac{F_{b,2}}{F_{b,1}}\right)^{-n_{2}}\left(\frac{F_{b,3}}{F_{b,2}}\right)^{-n_{3}}\left(\frac{F}{F_{b,3}}\right)^{-n_{4}},&F_{b,3}>F\geq F_{b,4}\\ \\ \ldots&\ldots\\ \\ \left[\prod_{i=1}^{k-1}\left(\frac{F_{b,i+1}}{F_{b,i}}\right)^{-n_{i+1}}\right]\left(\frac{F}{F_{b,k}}\right)^{-n_{k+1}},&F_{b,k}>F\end{array}\right..$$

$$\frac{d{N}_p}{dS}(S;\mathbf{\theta })=\frac{1}{E_p}\frac{d{N}_p}{dF}(F=S/E_p;\mathbf{\theta }).$$

$$P_k^{(p)}=\frac{1}{k!}{\frac{d^k{\mathcal{P}}^{(p})(t)}{d{t}^k}|}_{t=0}.$$

$${\mathcal{P}}_{\mathrm{NP}}(t;\mathbf{\theta })=\prod _p\exp \left[\sum _{m=1}^{\infty }{x}_{p,m}(\mathbf{\theta })(t^m-1)\right],$$

$$x_{p,m}(\mathbf{\theta })={\int }_0^{\infty }dS\frac{d{N}_p}{dS}(S;\mathbf{\theta }){\int }_0^{1}df\rho (f)\frac{(fS{)}^m}{m!}{e}^{-fS}.$$

$${dN_{p}(S)\over dS}=A_{p}\left\{\begin{array}[]{cc}\left({S\over S_{b}}\right)^{-n_{1}}&S\geq S_{b}\\ \left({S\over S_{b}}\right)^{-n_{2}}&S<S_{b}\,,\end{array}\right.$$

$$\frac{dN_{p}}{dS}=A_{p}\left\{\begin{array}[]{lc}\left(\frac{S}{S_{b_{1}}}\right)^{-n_{1}},&S\geq S_{b_{1}}\\ \left(\frac{S}{S_{b_{1}}}\right)^{-n_{2}},&S_{b_{1}}>S\geq S_{b_{2}}\\ \left(\frac{S_{b_{2}}}{S_{b_{1}}}\right)^{-n_{2}}\left(\frac{S}{S_{b_{2}}}\right)^{-n_{3}},&S_{b_{2}}>S\end{array}\right.$$

$$\rho (r)={\rho }_0\frac{(r/r_s{)}^{-\gamma }}{(1+r/r_s{)}^{3-\gamma }},$$

$$J(\varphi )={\int }_{l.o.s.}{\rho }^2(r)ds$$

$$n(z,R)\propto \exp \left[\frac{-R}{5\text{leavevmode}\text{kpc}}\right]\exp \left[\frac{-|z|}{1\text{leavevmode}\text{kpc}}\right],$$

$${ϵ}_p\equiv 1-\text{CDF}[{\mu }_p,n_p],$$

Full text

††preprint: MIT-CTP/5104

Introduction. There has been an extensive debate in the literature over the origins of the Galactic Center Excess (GCE), an extended and roughly spherically symmetric gamma-ray source filling the region within $\sim 1.5$ kpc of the Galactic Center (GC), with energy spectrum peaking at $1-3$ GeV Goodenough and Hooper (2009); Hooper and Goodenough (2011); Hooper and Linden (2011); Hooper and Slatyer (2013); Daylan et al. (2016); Calore et al. (2015); Ajello et al. (2016). The leading hypotheses are a new population of unresolved gamma-ray pulsars, individually too faint to be detected but in aggregate yielding the excess Abazajian et al. (2014); Abazajian and Kaplinghat (2012); Hooper et al. (2013); Mirabal (2013); Calore et al. (2014); Cholis et al. (2015); Yuan and Ioka (2015); O’Leary et al. (2015); Ploeg et al. (2017); Hooper and Linden (2018); Bartels et al. (2018a, b), or alternatively a signal from annihilating dark matter (DM) (e.g. Goodenough and Hooper (2009); Daylan et al. (2016); Karwin et al. (2017)). The latter explanation, if confirmed, would be of extraordinary importance for our understanding of the Universe, as the first non-gravitational probe of the properties of DM and its interactions with visible particles.

The hypothesis that unresolved point sources (PSs) generate much or all of the GCE can potentially be distinguished from the DM hypothesis via photon statistics Malyshev and Hogg (2011); Lee et al. (2015). The DM signal is expected to be dominated by the smooth Galactic halo (although see Agrawal and Randall (2017)), and the probability of seeing a certain number of photons from a spatial pixel is obtained from the Poisson distribution based on the expected number of photons. For a population of unresolved sources, the positions of the individual sources are not known, and the scatter of the number of sources (of a given brightness) in a pixel must be taken into account. Qualitatively, a population of unresolved PSs will generally have a greater probability to generate pixels with a large number of counts (due to the presence of one or more sources) or a very small number of counts (due to an absence of sources) compared to an extended diffuse source with the same overall expected number of photons.

Standard “template fitting” methods, where the sky is modeled as a linear combination of components with distinct spatial morphologies, can be adapted to incorporate these differences in statistical behavior. Templates for diffuse components consist of an expected number of photons in each pixel, possibly with a free overall normalization factor; templates for populations of unresolved PSs are characterized additionally by the “source count function” (SCF), which describes the probability that a given source has a certain brightness (i.e. produces a certain expected number of photons). It is then possible to calculate the probability to observe a certain number of photons in each pixel, as a function of the coefficients of the various templates and the source-count function parameters, and to study the resulting overall likelihood as a function of these parameters. This approach is called non-Poissonian template fitting (NPTF) Malyshev and Hogg (2011); Lee et al. (2015, 2016).

The NPTF method has previously been used to study the inner Galaxy in gamma rays, modeling the sky as a linear combination of Galactic diffuse emission, the large structures known as the Fermi Bubbles, isotropic extragalactic diffuse emission, isotropically distributed extragalactic PSs, Galactic PSs tracing the disk of the Milky Way, and the GCE. Modeling the GCE as a linear combination of a DM signal and an unresolved point-source population, with identical spatial morphologies for the signal but differing statistics, Ref. Lee et al. (2016) found that there was a strong statistical preference for the presence of the point-source population, and that the DM contribution was consistent with zero.

In this Letter, we explore the robustness of the NPTF to the presence of additional physical contributions to the gamma-ray data, which are not captured by the standard choice of templates. Some such modifications – such as changing the exact morphology of the disk-correlated or GCE-correlated PSs – have already been explored Lee et al. (2016) and shown not to qualitatively change the preference for PSs in the GCE. Motivated by the results of wavelet studies that find enhanced small-scale power in this region Bartels et al. (2016); Balaji et al. (2018), we focus on the possibility that there could be Galactic PSs present that are not physically associated with the Galactic disk or the GCE. If such small-scale power is present but corresponds to a source population unrelated to the GCE, we hypothesize that the template-based techniques employed in Ref. Lee et al. (2016) could incorrectly attribute these PSs to the GCE, and consequently bias the reconstruction of any DM contribution.

We first explore this possibility by simulating a proof-of-principle scenario, where one set of templates is used to construct the simulated data, and these data are then fitted with a different set of templates, similar to the NPTF pipeline in Ref. Lee et al. (2016). This allows us to directly observe and explore a bias reducing the reconstructed DM contribution, in a case where all inputs are known. In particular, we show that this bias can occur, can drive the reconstructed DM contribution to a value consistent with zero and inconsistent with its true simulated value, and that this remains true even when an additional DM signal is injected into the simulated data. We then examine real Fermi gamma-ray data, where the distribution of the underlying source populations is not known. We demonstrate that when an additional simulated DM signal is injected into the real data, it is also (incorrectly) attributed to PSs, qualitatively similar to the behavior we expect in the biased case from our proof-of-principle example.

Methodology and Data Selection. To employ the NPTF method, we use the NPTF package NPTFit Mishra-Sharma et al. (2017), interfaced with the Bayesian interference tool MultiNest Feroz et al. (2009). The total number of live points for all MultiNest runs is nlive = 500. Simulated contributions from Poissonian templates are generated by taking a random Poisson draw in Python of the product of the template and a predetermined normalization. Mock data for PS populations are generated using NPTFit-Sim Rodd and Toomey . We analyze eight years of Fermi Pass 8 gamma-ray data, for energies $2-20$ GeV, restricting to the top quartile of events graded by angular reconstruction. Our region of interest (ROI) is within a $30^{\circ}$ radius of the GC, with a $|b|<2^{\circ}$ plane mask, where $b$ is Galactic latitude. PSs contained in the 3FGL catalog Acero et al. (2015) are masked at a $99\%$ C.L. containment radius of 0.778 degrees for 2 GeV photons. Further details are supplied in the supplementary material.

We define the “standard NPTF pipeline” to be a fit including Poissonian templates for the Galactic diffuse emission (“Diffuse”), isotropic emission (“Iso”), the Fermi Bubbles (“Bub”), and the GCE (“DM”), and non-Poissonian templates for PSs tracing the Galactic disk (“Disk PS”), isotropic emission (“Iso PS”), and the GCE (“NFW PS”); template definitions are included in the supplementary material. “NFW” here stands for Navarro-Frenk-White Navarro et al. (1996), a commonly-employed prescription for the DM density profile that has previously been used in modeling the GCE. We perform a fit on the real data using this standard pipeline; we calculate the medians of the posterior distributions for the various model parameters, and simulate data based on those parameter values. We also define a “baseline NPTF pipeline” which contains the same templates except that it omits the NFW PS template; the Bayes factor between the standard and baseline fits thus describes the strength of evidence in favor of GCE-distributed PSs.

Proof-of-Principle Example of the Impact of Unmodeled Point Sources. We first investigate the performance of the NPTF pipeline on simulated data, in the case where the data contain contributions that are not well-modeled by the standard templates (in the supplementary material we explore the case where the mock data matches the templates). We will show that it is possible to induce significant biases to the reconstructed parameters, depending on the additional contributions. As an example, we consider unmodeled populations of PSs tracing the Fermi Bubbles. Small gas clumps have recently been identified Di Teodoro et al. (2018) in the Bubbles, and this could generate apparent gamma-ray PSs if the gas density in the clumps is sufficiently high.

More generally, the Bubbles serve as an example of a component with a spatial morphology that is not degenerate with the disk or with isotropic PSs, and which could conceivably be absorbed by a GCE template in preference to the other PS templates. As we show in separate work Leane and Slatyer , we do not detect PSs associated with this Bubbles template in real data. As such, we do not expect that this example precisely describes the real gamma-ray sky; it is simply a proof-of-principle example that (as we will see) the reconstructed DM-associated flux can be significantly biased if the templates used in the fit do not adequately describe the real data.

Studying the Fermi data with the standard pipeline yields simulated data with essentially no DM contribution, consistent with the result of Ref. Lee et al. (2016). In order to simulate a DM contribution of comparable size to the GCE, we also perform a fit with the NFW PS template removed from the standard pipeline, and take the posterior median of the DM template normalization in this case. This gives a DM flux of about 2 percent of the total photon flux in the masked ROI. To simulate a hypothetical source population tracing the Fermi Bubbles, we use the same source count function parameters as simulated for the NFW PS population.

We simulate PSs correlated with the Bubbles and Disk, as well as a DM component and diffuse Bubbles, Isotropic and Galactic Diffuse templates. We then apply the standard NPTF pipeline to this scenario (minus the omitted isotropic PS template); in particular, the fit includes a template for NFW PSs, but not Bubbles-correlated PSs. We also repeat the same fit with the analogue of the baseline NPTF pipeline for this case (again omitting the isotropic PS template).

Table 1 summarizes the relative Bayes factors and flux containment intervals in our analyses. Comparing the Bayesian evidences for these two models yields the Bayes factor in favor of NFW PSs; we consistently find a preference in favor of the model containing NFW PSs, despite the fact that no such PSs are present in the simulation.

Figure 1 shows the best-fit fluxes attributed to the various templates in the “standard pipeline” case, as a fraction of the total simulated photon flux in the ROI. These are the first of our main results: The fit misattributes the simulated DM signal to point sources. The DM-sourced photon flux appears to be divided between NFW PSs and smooth emission in the Bubbles. Increasing the amount of injected DM (but keeping the other simulated components the same), the DM flux continues to be absorbed by the NFW PS template, until it reaches a threshold of absorption, where some DM is eventually identified.

We note that a bias-inducing population of unmodeled PSs could also be consistent with wavelet studies indicating an excess of small-scale power at the GCE, exceeding expectations from the Galactic diffuse emission; we leave a detailed study of consistency to future work. Non-NPTF studies Bartels et al. (2016) have argued that this small-scale power cannot be easily explained by the unresolved PSs associated with known populations, but have not been able to clearly associate it spatially with the GCE. Consequently, their results also do not place a firm upper limit on any DM contribution to the GCE.

Dark Matter Injection Test in Real Gamma-Ray Data. Above we have identified a proof-of-principle example of how a DM signal can be misattributed to PS by the NPTF, when an unmodeled PS population is present. We now seek to test whether similar effects could be occurring in the real gamma-ray data.

The complicating factor is that in real data, the distributions of any unmodeled PS populations are (by definition) unknown. One possible test is to add additional physically-motivated templates for PS populations, to see if the GCE is then attributed to a DM component, but we have not yet identified an additional template that has this effect (and as the number of templates tested proliferates, caution would be needed in statistical interpretation of the results).

However, as a possible diagnostic that does not rely on specifying the biasing population, we can artificially inject a DM contribution into the data and see if it is correctly reconstructed. In our proof-of-principle example, increasing the DM contribution well above the GCE level still resulted in the attribution of the DM component to NFW PSs. It would be a coincidence if, in the presence of an unmodeled PS population (or another effect that biased the normalization of the DM component), the bias was sufficient to incorrectly attribute the observed GCE to PSs, but could not absorb any additional DM injections. In other words, if this effect is present, there is no reason to think it should be saturated.

As previously, we simulate this additional DM injection by taking a random Poisson draw in Python of the DM template multiplied by an appropriate normalization factor; the normalization is varied to test a range of injection levels.

Table 2 presents our results for injected DM fluxes as a fraction of the total photon flux in the real + simulated data in the ROI (post-injection); the reconstructed fluxes associated with various templates, within the standard and baseline NPTF pipelines; and the Bayes factors between pairs of models. We also perform fits where the DM template is fixed to the injected value, to test the Bayes factor for this scenario compared to the cases where DM is allowed to float freely.

In general, we find that when the injected DM signal is up to a factor of $\sim 5$ larger than the baseline GCE, the injected signal is completely attributed to the NFW PSs template, with the reconstructed bounds on the DM signal being consistent with zero flux and strongly inconsistent with the injected value. We find that the model where the NFW PS and DM templates are allowed to float, with the DM template being driven toward zero flux, is preferred over both the model where the DM contribution is fixed to its injected value, and (to a larger degree) the model without NFW PSs. This behavior – where a model that is inconsistent with the injected signal is favored by the fit – is qualitatively similar to the simulated proof-of-principle case, and consistent with a mismodeling effect capable of driving any actual DM contribution to the GCE to zero flux.

Once the injected signal is sufficiently large, the flux associated with the DM template becomes inconsistent with zero, albeit still inconsistent with the injected value; however, this requires a very large injection. For example, injecting a DM signal almost an order of magnitude brighter (15.2$\%$ post-injection) than the GCE itself results in around half the injection (only $\sim 8\%$) being attributed to DM.

Figure 2 shows results for recovered flux fractions in the ROI for varying levels of DM injection: 0, 1.8, 6.7, and 15.2$\%$ (post-injection). The top-left plot shows results of a fit on the real data, without any injected DM signal, to serve as a baseline for comparison. By comparison with the no-injection case, we see two important effects: firstly that (as noted above) the flux attributed to the DM template is consistent with zero and inconsistent with the injected value, and secondly that the NFW PS flux fraction increases, approximately absorbing the injected DM signal. As the DM injection amount increases, we see that the NFW PS flux fraction continues to increase, until it reaches a saturation point and the DM template begins to absorb some of the flux. In order for the DM to be detected with non-zero flux, the injected DM signal appears to require a total flux a factor $\gtrsim 5$ larger than the GCE itself.

Conclusions and Outlook. We have studied examples of how NPTF methods can be biased in both real and simulated gamma-ray data, and how this could impact explanations of the GCE. We have showed a proof-of-principle example in simulated data where a DM signal can incorrectly be attributed to PSs by the NPTF, as a result of PSs with a spatial distribution that is not described by the standard templates.

We have searched for behavior consistent with this effect in the real Fermi gamma-ray data. If an unmodeled source population is biasing results of NPTF studies, causing the NFW PS template to absorb the DM flux, it is likely not saturated in its ability to do so. We injected an artificial DM signal into the real data, with various normalizations, to test if similar behavior is observed to the simulated-data proof-of-principle example. Strikingly, the injected DM signal was misattributed to PS, with the recovered NFW PS flux fraction increasing proportionally to the size of the DM signal injected. We increased the normalization of the DM signal, and found this effect persists until the amount of DM injected exceeds the GCE itself by almost an order of magnitude, at which point the fit recovers a DM component, but not with the correct injected normalization. This behavior is qualitatively similar to the simulated proof-of-principle example; in both cases there are strong statistical preferences for models that reconstruct too little dark matter, compared to the injected values. This suggests that the data are mis-modeled in some way, and that DM may be the dominant contribution to the GCE after all - although this does not constitute positive evidence for DM, and the GCE could still be composed mainly of PSs. In the supplementary materials, we briefly explore the degree to which this evidence of mis-modeling persists across different diffuse models and ROIs.

An improved future version of the NPTF pipeline, with greater flexibility to accommodate a diverse range of source populations or a different range of ROIs, could potentially characterize and eliminate biases similar to those demonstrated in our examples. Upcoming observations by radio telescopes such as MeerKat and the Square Kilometer Array (SKA) should allow some currently unresolved source populations to be identified Calore et al. (2016). Complementary studies such as those using wavelet techniques Bartels et al. (2016); Balaji et al. (2018), or focusing on improved diffuse emission modeling with new frameworks such as SkyFact Storm et al. (2017), could also help disentangle sources of emission from the GC.

Acknowledgments. We thank Tim Linden, Mariangela Lisanti, Nick Rodd, and Ben Safdi for helpful discussions and comments. We thank Tim Linden for suggesting the injection test on real data, the authors of Ref. Macias et al. (2018) for providing their Boxy Bulge template, and the Fermi Collaboration for the use of Fermi public data. RKL thanks the DESY Theory Group and the Princeton Institute for Advanced Study for their hospitality during this work. RKL and TRS are supported by the Office of High Energy Physics of the U.S. Department of Energy under Grant No. DE-SC00012567 and DE-SC0013999; TRS is partially supported by a John N. Bahcall Fellowship.

I Template Fitting Methodology

In this section we review the template fitting methods employed in this work, describe the templates we use in our analysis, and detail the regions of interest and dataset studied.

I.1 Poissonian Template Fitting

Template fitting is a widely used approach that fits a linear combination of templates for physical sources of photon emission (in this case, gamma rays) to observational data. In general, templates may vary with both position and energy, but for our purposes we will work within a single broad energy bin and characterize templates purely by their spatial distributions (and photon statistics).

In the simplest scenario, such spatial templates describe the expected mean photon counts per pixel, photons are assumed to be independent events, and thus the probability of observing a certain number of photons sourced by a given template in a specific pixel is given by the Poisson distribution. This assumption of independence is valid for diffuse emission and for PS populations where the locations of the PSs are known.

Contributions from different templates can be trivially combined in this case. Let us denote the pixel-dependent expected photon counts from the individual spatial templates by $\mu_{p,\ell}({\bf\theta})$, where $p$ is a pixel index, $\ell$ indexes the choice of template, and ${\bf\theta}$ is a set of parameters controlling the properties of the templates. In the Poissonian case, we will generally choose the parameters ${\bf\theta}$ to control the coefficients of the various templates, so that:

[TABLE]

where $T_{p,\ell}$ defines a spatial template with some fixed normalization. The mean number of total counts within a specified pixel, $\mu_{p}({\bm{\theta}})$, is then given by:

[TABLE]

and the probability of drawing $n_{p}$ counts in a given pixel $p$, as a function of the parameters $\bm{\theta}$, is:

[TABLE]

In this case, for given model parameters $\bm{\theta}$, a model $\mathcal{M}$ can be fit to a dataset $d$ with $n_{p}$ counts in each pixel $p$. The likelihood function is then just the product of Poisson likelihoods,

[TABLE]

I.2 Non-Poissonian Template Fitting

In the situation where the gamma-ray emission contains spatial correlations between photons that are not captured by the spatial template itself, the photons from a given pixel can no longer be treated as independent events, and Poissonian statistics are not appropriate. This is the case, for example, where there are PSs producing two or more photons in the data, but the positions of these sources are not included in the spatial template; they are unresolved. The presence of the sources manifests itself in an increased probability of seeing multiple photons from the same location, compared to Poissonian expectations; we can search for this effect while remaining agnostic as to the exact positions of the sources.

The detailed statistical treatment of this situation was first discussed in the context of gamma-ray data in Ref. Malyshev and Hogg (2011), and expanded in Ref. Lee et al. (2015) (other techniques for estimating contributions from unresolved PS populations are discussed in e.g. Zechlin et al. (2016a, b)). Ref. Lee et al. (2016) realized that these techniques could be combined with template-fitting methods to include information on the spatial distributions of the backgrounds and potential point-source populations, and tested the resulting NPTF pipeline on Fermi data; follow-up studies were conducted in Refs. Linden et al. (2016); Lisanti et al. (2016), and a public code package (which we employ) was presented in Ref. Mishra-Sharma et al. (2017). We summarize the key features of the mathematical framework below.

In the NPTF, the photon count probability distribution $P_{n_{p}}^{(p)}$ discussed above has an additional dependence on a pixel-dependent PS source-count distribution $dN_{p}/dF$, where $F$ is the source flux. We model the source count function as a multiply broken power law, and assume it is the same in all pixels up to a pixel-dependent normalization factor:

[TABLE]

Here, the source-count distribution is parameterized with an arbitrary number of breaks $k$, denoted by $F_{b,i}$ with $i\in[1,2,\ldots,k]$, and $k+1$ indices $n_{i}$ with $i\in[1,2,\ldots,k+1]$. The spatial distribution of the PSs is described by the template $T_{p}^{({\rm PS})}$, with a pixel-independent normalization $A({\bm{\theta}})$ which is a function of the model parameters $\bm{\theta}$. Note that while the number of sources varies between pixels, in our work the locations of the flux breaks and the indices are fixed across pixels.

In practice we need the source count function as a function of photon counts $S$, rather than photon flux $F$; the two are related by the pixel-dependent Fermi exposure map $E_{p}$. The conversion formula is:

[TABLE]

To calculate the pixel likelihoods for models containing templates with non-Poissonian templates, the method of generating functions is employed; the generating function $\mathcal{P}^{(p)}(t)$ for the complete model (consisting of a sum over templates) is the product of the individual template generating functions, and the probability $P_{k}^{(p)}$ to observe $k$ photons in pixel $p$ can be derived from this generating function as:

[TABLE]

Here $t$ is an auxiliary variable. The generating function for a Poissonian template is given by $\mathcal{P}_{\ell}^{(p)}(t)=e^{\mu_{p,\ell}(t-1)}$, while for the non-Poissonian PS templates we have instead:

[TABLE]

where

[TABLE]

The point spread function (PSF) of the Fermi instrument is taken into account through $\rho(f)$, which describes the dilution of the flux from a PS into neighboring pixels due to smearing by the PSF. For a distribution $dN_{p}(S;{\bm{\theta}})/dS$, $x_{p,m}$ describes the expected number of $m$-count PSs in the pixel $p$. For more details, see Refs. Malyshev and Hogg (2011); Lee et al. (2015).

We generate results using the public Python and Cython NPTF package NPTFit Mishra-Sharma et al. (2017), interfaced with the Bayesian interference tool MultiNest Feroz et al. (2009), which implements the likelihood framework described in this section. The total number of live points for all MultiNest runs is nlive = 500. Typically about $10^{6}$ posterior samples are generated for each run.

I.3 Data Selection

We use the Pass 8 Fermi data, which were collected in the energy range $2-20$ GeV over 413 weeks, from August 4th 2008 to July 7th 2016. The event class UltraCleanVeto (1024) was applied, which has the highest cosmic ray rejection. This is further restricted to the top quartile of events graded by angular reconstruction (PSF3 (32)), with quality cuts DATA_QUAL==1 && LAT_CONFIG==1. The maximum zenith angle is $90^{\circ}$.

I.4 Modeling and Masking of Resolved Sources

The gamma-ray sky contains a number of high-significance PSs. To mask and/or model these sources, we use as a reference the Fermi-LAT 3FGL catalog Acero et al. (2015). While there are more up-to-date catalogs, in particular the 1FIG catalog Ajello et al. (2016), the main goal of masking PSs is to eliminate very bright pixels where computing the relevant photon count probabilities is time-consuming, and the 3FGL catalog has proved adequate for this purpose (it has also been demonstrated that masking the 1FIG PSs does not affect the GCE Bartels et al. (2018c)).

I.5 Regions of Interest

We study two different regions of interest (ROI) in the inner Galaxy, one with known PSs masked, and one with no masking of PSs.

In detail, the regions are:

• •

Inner Galaxy (masked): within $30^{\circ}$ of the GC, and $|b|>2^{\circ}$, with all 3FGL PSs masked at a $99\%$ containment at an energy of 2 GeV (corresponding to a radius of 0.778 degrees), and

• •

Inner Galaxy (unmasked): within $30^{\circ}$ of the GC, and $|b|>2^{\circ}$,

where $\ell$ is longitude, and $b$ is latitude, in the galactic coordinate system. The unmasked case provides more ROI close to the GC, as the 3FGL PSs are clustered around the GC, and masking the 3FGL sources removes a substantial fraction of the near-GC solid angle. We perform all tests in the $30^{\circ}$-radius ROI with the 3FGL sources masked, and then check the key results with the unmasked ROI.

In principle, Fermi’s exposure should be accounted for individually in each pixel, when converting between the flux of a source and its expected counts as in Eq. S6. However, in practice this is not necessary, and we instead subdivide the ROI into a number of “exposure regions”. Within each such subregion, the exposure is approximated as its average value over the exposure region, for the purpose of converting between flux and counts. For the Poissonian templates, where the exposure simply modifies the expected number of counts in each pixel, we do not perform this approximation and instead implement the full exposure map. In our Inner Galaxy analyses, we break the sky into 5 exposure regions. To ensure this is an accurate description of the data, we have checked stability of our main results (discussed in the main text) up to at least 30 exposure regions.

I.6 Spatial Templates

By default we use the templates provided publicly with the NPTFit package. We review the properties of these templates below, and note when we differ from the NPTFit baseline.

All templates are implemented in a HEALPix Gorski et al. (2005) pixelization of the data with nside = 128, such that the photon map contains 196,608 equal area pixels over the full sky. All Poissonian templates are exposure-corrected; they are thus maps of expected photon counts/pixel (in contrast to non-Poissonian templates which are maps of sources/pixel). The Galactic diffuse emission model, which accounts for the bulk of the gamma-ray emission, is smoothed by the correct Fermi-LAT PSF using the Fermi Science Tools. The DM and Fermi Bubbles diffuse templates are smoothed with a Gaussian with $\sigma=0.1812$ degrees. All templates are then rescaled to have a mean of one in the region with $|b|>2^{\circ}$ within $30^{\circ}$ of the GC.

Table S1 details the priors used for each template parameter. When floating any template, we scan with a log-flat prior distribution for the template normalization. The remaining parameters have linear-flat prior distributions. We describe each spatial template below.

• •

Galactic Diffuse Emission

Diffuse gamma-ray emission from the Milky Way dominates the gamma-ray sky. Such emission originates from three main sources: (1) cosmic-ray (CR) protons colliding with the gas and producing photons via pion production, $pp\rightarrow X+\pi^{0}\rightarrow X+\gamma\gamma$, (2) CR electrons upscattering photons of the interstellar radiation field to gamma-ray energies via inverse Compton scattering, and (3) CR electrons scattering on the gas to produce photons via bremsstrahlung (this third process is generally subdominant to the other two).

To model this diffuse background emission, we use the Pass 6 Fermi diffuse model p6v11; this is the most recent diffuse model released by the Fermi Collaboration that does not incorporate a built-in template for the Fermi Bubbles or add a component for otherwise-unmodeled diffuse emission, allowing us to study these contributions separately. This was also the principal diffuse background model used in earlier NPTF studies of the GCE, facilitating comparisons Lee et al. (2016).

However, this diffuse model was fitted to much earlier data from Fermi, with different instrument response functions to the current dataset; to translate the flux model into counts we use the appropriate instrument response functions for our current dataset, but errors in the instrument response modeling (either now or in the construction of the p6v11 model) could in principle introduce difficult-to-model systematic uncertainties.

As a cross-check of our results, we also employ two models for the Galactic diffuse emission generated using GALPROP Strong and Moskalenko (1998), labeled models A and F in Ref. Calore et al. (2015). These models were found to provide a better description of the data than the p6v11 model at GeV+ energies. They contain two component templates which are floated independently, corresponding on one hand to gas-correlated emission from pion production and bremsstrahlung, and on the other to emission from inverse Compton scattering.

The Galactic Diffuse emission template is specified solely by its normalization value $A_{\rm dif}$, or by two normalization values $A_{\pi^{0}\text{brem}}$ and $A_{\rm ics}$ in the case of the GALPROP-based diffuse models.

• •

Fermi Bubbles

The Fermi Bubbles are a large gamma-ray structure extending nearly $\sim 10$ kpc on either side of the galactic disk Su et al. (2010), and approximately centered at the Galactic Center.

We use the Fermi Bubbles template from Ref. Su et al. (2010). This is taken to have uniform emission intensity before exposure correction. The smooth Bubbles template is characterized by only the normalization $A_{\rm bub}$, with log-flat prior distribution. We leave an analysis of possible alternate models of the Bubbles Ackermann et al. (2014) to later work.

Hypothetical PS emission from the Fermi Bubbles region is specified by extra model parameters, ${\bm{\theta}}=\{A_{\rm bub}^{\rm PS},n_{1}^{\rm bub},n_{2}^{\rm bub},S_{b}^{\rm bub}\}$. The PS source-count function in a given pixel $p$ is taken to be a singly broken power law:

[TABLE]

where $S_{b}=S_{b}^{\rm bub}$ is the break, $n_{1}=n_{1}^{\rm bub}$ and $n_{2}=n_{2}^{\rm bub}$ are the slopes above and below the break, and $A_{p}=A_{p}^{\rm bub}$ is the pixel-dependent normalization. In order for the total number of photons from PS to be finite, it is required that $n_{1}$ and $n_{2}$ are greater than and less than 2 respectively. The source-count function in counts in Eq. (S10) is related to the function in flux via Eq. (S6).

• •

Isotropic Background

Gamma-ray emission from extragalactic sources is expected to be approximately isotropic on the sky; this component can also absorb residual cosmic-ray contamination (although in the UltraCleanVeto dataset we use, such contamination should be very small). We include both a Poissonian isotropic component, specified by its normalization $A_{\rm iso}$ (with log-flat prior), and also a component for isotropically distributed PSs. In high-latitude studies of the isotropic PS population, it has been shown that the choice of priors can determine the source-count function at low flux Lisanti et al. (2016), and that more than one break is required to accurately describe the source-count function at high flux. Consequently, we model the isotropic PS source count function with two breaks, so it is specified by the model parameters ${\bm{\theta}}=\{A_{\rm iso}^{\rm PS},n_{1}^{\rm iso},n_{2}^{\rm iso},n_{3}^{\rm iso},S_{b_{1}}^{\rm iso},S_{b_{2}}^{\rm iso}\}$.

Explicitly, the PS source-count function in a given pixel $p$ is a doubly broken power law:

[TABLE]

where $S_{b_{1}}=S_{b_{1}}^{\rm iso}$ is the first break, $S_{b_{2}}=S_{b_{2}}^{\rm iso}$ is the second break, $n_{1}$ and $n_{2}$ are the slopes above and below the first break, $n_{2}$ and $n_{3}$ are the slopes above and below the second break, and $A_{p}=A_{p}^{\rm iso}$ is the pixel-dependent normalization. The NP model priors other than normalization are taken to be linear flat. The source-count function in counts in Eq. (S11) is related to the function in flux via Eq. (S6).

• •

Generalized NFW Profile (Galactic Center Excess)

The flux profile of the GCE has been previously found to be well described by the square of a generalized NFW profile Daylan et al. (2016); Calore et al. (2015), consistent with DM annihilation if the DM density profile follows the (generalized) NFW profile. As such, we include a template for a signal with spatial distribution matching that of a squared generalized NFW profile Navarro et al. (1996) integrated along the line of sight (l.o.s.).

The functional form of the density profile is:

[TABLE]

where $r_{s}=20$ kpc, $\gamma=1.25$, and the normalization is allowed to float. (We also check the impact of including a template of a standard NFW profile with $\gamma=1.0$, and find no significant difference in our results.) The resulting flux is determined by:

[TABLE]

where $s$ is the l.o.s. distance and $\phi$ is the angle from the GC.

The Poissonian generalized-NFW template, as would be appropriate for a DM signal, is characterized by its normalization $A_{\rm NFW}$. PSs distributed following this profile are characterized by the model parameters ${\bm{\theta}}=\{A_{\rm NFW}^{\rm PS},n_{1}^{\rm NFW},n_{2}^{\rm NFW},S_{b}^{\rm NFW}\}$, which determine the source count function via Eq. (S10). We will refer to “DM” to mean the generalized NFW profile ($\gamma=1.25$) template, rather than the traditional NFW profile ($\gamma=1.0$).

• •

Galactic Disk

We model PS emission from the galactic disk through a doubly exponential thin-disk source distribution,

[TABLE]

with radius $R=5$ kpc and a scale height z=1\leavevmode\nobreak\kpc. (This disk model, and a thinner disk with scale height z=0.3\leavevmode\nobreak\kpc, were both tested in Ref. Lee et al. (2016); the resulting differences in the fit were found to be small.) We integrate along the line of sight to obtain the expected number of sources per pixel.

PSs distributed along the disk are characterized by the model parameters ${\bm{\theta}}=\{A_{\rm disk}^{\rm PS},n_{1}^{\rm disk},n_{2}^{\rm disk},S_{b}^{\rm disk}\}$, which determine the source count function via Eq. (S10).

• •

Boxy Bulge

The Boxy Bulge is an X-shaped structure at the center of the Milky Way, which appears to be correlated with gamma-ray emission from the inner Galaxy Bartels et al. (2018d); Macias et al. (2018, 2019); Song et al. (2019). We do not use the corresponding template in the standard or baseline pipeline, but since it has been suggested that stars associated with the Boxy Bulge could contribute significantly to the GCE, we will later perform the test of replacing the NFW PS template with a PS population tracing this structure.

We use the template for the Boxy Bulge from Ref. Macias et al. (2018), and pixelize the template with nside=128. Smooth emission from the Boxy Bulge template is characterized by the normalization parameter $A_{\rm bul}$. PS emission associated with the Boxy Bulge is described by the model parameters ${\bm{\theta}}=\{A_{\rm bul}^{\rm PS},n_{1}^{\rm bul},n_{2}^{\rm bul},S_{b}^{\rm bul}\}$, which determine the source count function via Eq. (S10).

II Proof-of-Principle with Simulated Data: Supplementary Figures and Analyses

In this section, we supplement the main text by checking that we correctly recover the simulated templates, when the same templates are used in the construction of simulated data and at the analysis stage.

We also study simulated data biases in a ROI without a PS mask (in contrast to the masked ROI in the main text) for the example of a PS population tracing the Fermi Bubbles. In this case, we also include isotropic PSs in the fit and the resulting simulated data. In the previous (masked) analyses, the masking of the 3FGL PSs ensured that the only isotropic PSs in the data were faint and unresolved, so the effect of leaving them out was expected to be relatively minor; in contrast, in the unmasked case, the isotropic PS population is expected to include some number of very bright sources.

Table S2 details all the model parameters used to create the simulated data, for the two inner Galaxy ROIs (with 3FGL PSs masked and unmasked, respectively). Note that not all templates are simulated in every analysis.

II.1 Recovery of Parameters from Simulated Data

We simulate mock data containing disk PSs, bubbles PSs, DM, and all other Poissonian templates contained in the standard pipeline (including isotropic PS in the unmasked ROI). We test that analyzing this simulated data with the same pipeline correctly recovers the simulated parameters within the uncertainties (for this realization), and in particular correctly identifies the smooth DM component as well as the two PS templates.

Figure S1 shows the source count functions and flux fractions attributed to the various templates when the templates used in the simulation match those used in the analysis. We also check if increasing the normalization of the simulated DM component impacts the analyzed results. We find it correctly recovers all simulated parameter values and fluxes in all cases.

II.2 Supplementary Results for Simulated Point Sources in the Fermi Bubbles

In addition to the results in the main text, corresponding to the ($30^{\circ}$ radius, $|b|<2^{\circ}$ plane mask) masked ROI, we test the same proof-of-principle example in the ($30^{\circ}$ radius, $|b|<2^{\circ}$ plane mask) unmasked ROI.

Table S3 shows the result of applying the standard and baseline NPTF pipelines to simulated data containing DM and Bubbles PSs in this case, when an additional DM component is injected. As in the masked ROI, we find that the fit incorrectly reconstructs the simulated DM component as NFW PSs when given the option.

Figure S2 shows the posterior probability distributions for the relevant template flux fractions, and the reconstructed source count functions for the PSs.

Figure S3 shows, as a supplement to the main text, the evolution in the reconstructed source count functions for disk and NFW PSs, for 4 different injected DM fluxes (2, 5, 7, and $15\%$ post-injection), in the masked ROI. These results correspond to the flux posteriors shown in Fig. 1, where the mock data contains Bubbles PSs and a DM component, but is fitted with the standard NPTF pipeline (i.e. with NFW PSs and without Bubbles PSs). We observe that as the DM normalization is increased, the source count function for the reconstructed NFW PSs remains fairly stable; there is some trend toward higher $dN/dF$ at low $F$, but the difference remains within the uncertainty bands.

Note that in general for our Bubbles analyses, as the Fermi Bubbles extend to much higher latitudes than our ROIs, it is also possible to fix the diffuse emission within the Bubbles to its high-latitude normalization, instead of allowing it to float. In this case, we find that even more flux is attributed to the NFW PSs. Similarly, fixing the diffuse model to its high-latitude value provides qualitatively comparable results.

We generate 20 sets of simulated data, and find this result persists, with comparable posteriors and Bayes factors.

It is possible that biases of this type could be detected and modeled by an improved version of the NPTF pipeline. For example, giving extra spatial freedom to the NFW PS template, or allowing the source count function to vary between different spatial subregions of the template, could identify a mismatch between the NFW template and the true PS distribution. Varying the ROI and comparing the results to those expected in simulated data could also help distinguish scenarios with different spatial PS distributions. For example, the Bayes factor in favor of NFW PSs may change between scenarios when the resolved 3FGL sources are unmasked or masked, and may not change in a manner matching simulations. Reducing the ROI to the region where the GCE is brightest can similarly help eliminate biases from populations outside that region. Both strategies, however – increasing the number of degrees of freedom, and reducing the size of the ROI – are also expected to inflate the statistical uncertainties, reducing constraining power. We have studied some possible improvements to the NPTF pipeline in this section and the following sections, but leave a comprehensive study to future work.

III Injection test on real data: Supplementary Figures and Analyses

Figure S4 shows, as a supplement to the main text, the evolution in the reconstructed source count functions for disk, isotropic and NFW PSs, for no injection as well as 3 different injected DM fluxes (1.8, 6.7, and $15.2\%$ post-injection), in the masked ROI. These results correspond to the flux posteriors shown in Fig. 2.

Table S4 shows the reconstructed DM flux fractions for different combinations of templates, and the relative Bayes factors for these scenarios, in the unmasked ROI, for a single normalization of the injected DM signal (corresponding to an injection of similar size to the GCE). As previously, we see that the DM component is reconstructed to a value consistent with zero and inconsistent with the injected normalization. The scenario where the DM normalization is floated is preferred over the case where it is fixed to the injected value, and also strongly preferred over the case with no NFW PSs. This behavior is again qualitatively consistent with the results we see in the proof-of-principle example using simulated data.

Figure S5 shows the posterior probability distributions for the template flux fractions in the case of a $1.8\%$ DM injection when the 3FGL are unmasked, as well as the source count functions for the reconstructed PS populations.

In Figure S6, we show for purposes of comparison the recovered flux fractions and intensities when the Fermi data, plus an injected DM signal (with normalization equal to the best-fit DM signal under the baseline pipeline), is analyzed with the baseline NPTF pipeline. In this case, there is no NFW PS template available to absorb the GCE, and both the original GCE and the injected signal are attributed to the DM template, which rises in normalization accordingly.

Note to ensure we are not just observing an extreme statistical fluctuation in the simulated DM signal injected into the real data, for the $1.8\%$ injection we also generate 100 sets of simulated DM signals, which are each independently redrawn from a Poisson distribution. We then analyze 100 independent datasets of the Fermi data injected with the independent simulated DM signals. We find that the recovered DM flux is always significantly below the true injected value, and is consistent with the values presented in Tab. 2. As such, the DM fluxes recovered and Bayes factors as shown are representative values. (Note the containment bands for Bayes factors for all analyses in Tab. 2, based on 20 sets of independent simulations and analyses, are shown in Fig. S8).

We also check the impact of injecting a DM signal into an analysis where NFW PSs are interchanged in the fit with Boxy Bulge PSs, to test the possibility that the behavior of the real data under DM injection could simply be due to a mismodeling of PSs associated with the Boxy Bulge. We find that the Boxy Bulge PS template behaves very similarly to the NFW PS template; in particular, it similarly absorbs the injected DM flux.

III.1 Dark Matter Injection into a Simulated Version of the Best-Fit to the Fermi Data

We also consider the case where the best fit values using the standard pipeline on the real data with no injection are used to produce simulated data. We then inject this simulated data with the exact same DM normalizations as Fig. 2, and then analyze this data with the exact same templates as simulated. This serves to highlight that the data do not appear to be accurately described by the standard templates – in this scenario fluxes are approximately recovered within uncertainties. Specifically, of the 20 realizations, we observe for $1.8\%$ flux injection that 7 out of 20 simulated realizations recover the correct flux within the $68\%$ containment bands, 14 of 20 for $6.7\%$ injection, and 14 out of 20 for $15.2\%$ injection. This is approximately consistent with expectations, although the reconstruction appears to be biased slightly low for the $1.8\%$ injection case. Importantly, in the case where the DM injection is comparable to the GCE itself, there is no statistical preference for the scenario where the DM contribution is allowed to float (and is reconstructed at a low value) over the scenario where the DM contribution is fixed to its injected value, unlike what is found with the injection test on real data.

In the case with no DM injection, the Bayes factor is similar to what we observe in our proof-of-principle example – although this similarity should be viewed with some caution, as for simplicity, our simulated data in the (masked) proof-of-principle case did not include the isotropic PS population. Nonetheless, this suggests that it could be difficult to distinguish between a case similar to our proof-of-principle scenario, where a real DM signal is incorrectly reconstructed as PSs due to the presence of an unmodeled PS population, and the case where the GCE is entirely PSs, solely via examination of the Bayes factor in favor of PSs. This contrasts with the case where the GCE receives a $50\%$ contribution from DM and $50\%$ from PSs, studied in Ref. Lee et al. (2016), where it was found that while the GCE could be incorrectly reconstructed as originating solely from PSs, the Bayes factor in favor of PSs was generically much smaller.

Figure S7 shows the recovered flux posteriors in this scenario, in one particular representative realization. Compared to Fig. 2, we observe that the injected DM and other simulated templates are much closer to being successfully recovered.

Table S5 details the reconstructed DM flux fractions for different combinations of templates, and the relative Bayes factors for these scenarios, in the same realization.

Figure S8 presents the difference in the spread of Bayes factors between the simulated data based on the best-fit model with the standard pipeline, and the real data with and without injection. This again highlights the data do not appear to be described accurately by the standard templates.

III.2 Increasing the Number of Breaks in the NFW Source-Count Function

Very faint PSs approach the limit of statistically behaving just like DM: their probability to emit $n$ photons in a given pixel becomes effectively zero for $n>1$, as the probability to emit observed photons at all is so low. As such, it is conceivable that the reason the DM flux is not being recovered is because it is being injected into the soft end of an overly-constrained singly-broken power law SCF – for example, the behavior of the SCF might be almost entirely dictated by the bright sources which are near the single break. As such, we test the effect of enforcing a steep slope right below the 1-photon threshold, as well as allowing a break around this threshold (prior range of the additional break taken as $S_{b}=[0.1,1]$), with a slope below this threshold that is allowed to float.

Figure S9 shows how these modifications change the recovered SCF for the case with 3FGL PSs masked and a 1.8% DM flux injected into the data (though we have also checked we recover the same behavior unmasked).

We find no evidence that such a break is preferred by the data; whether the slope below the break is floated or forced to be steep, we find the change in log evidence relative to the baseline case is within the error bars of the log-evidence, and the flux posteriors for the various components do not change.

IV Preference For the Dark Matter Normalization to Float Negative

Prior to this section, we have always forced the normalizations of all templates to be positive, since they represent physical emission mechanisms; we will retain this convention unless explicitly noted otherwise. However, as an alternative way to understand the behavior of the analysis when additional DM is injected, we can expand the prior on the coefficient of the DM template so that it includes negative values. This is clearly not a realistic scenario, but can indicate the presence of systematic effects which drive down the reconstructed DM component. Analyzing the real data (with no additional DM injection), with the p6v11 Galactic diffuse emission model, we find that the DM component prefers to float to a very negative value when given the option.

Figure S10 demonstrates this result; the DM component is reconstructed with a negative flux $\sim 4-5\times$ larger in magnitude than the baseline GCE. The NFW PS template rises in compensation, such that the sum of their fluxes remains approximately constant. This provides an intuitive understanding of the results of the injection test; adding an extra DM signal to the data reduces the degree to which the DM component prefers to run negative, but a positive signal will only be reconstructed once enough simulated DM has been added to completely cancel out the preferred negative flux.

Furthermore, until this threshold is reached, in fits where the DM component is forced to be non-negative, both the DM coefficient and the amount of power in very faint NFW PSs (degenerate with the DM component) will be forced toward zero. Due to the degeneracy between these components, only the sum of their fluxes is well-constrained, but if their sum prefers a negative value and each component is forced to be non-negative, the unique best-fit solution will set both contributions to zero. Consequently, we expect (as observed) both a near-zero reconstructed DM component and a source count function for the NFW PSs that has very little power at low flux. Adding additional freedom to the source-count function at low flux, as studied above, is not expected to change the results of the analysis, as in this case the degeneracy between DM and low-flux PSs has already been broken (by both components being forced to zero). Only once sufficient DM is injected (or when systematic factors causing the oversubtraction of the smooth component are corrected), and the total preferred power in DM and low-flux sources becomes positive, should the degeneracy between the two components become a potentially important systematic issue.

We can also ask whether this behavior, where the DM coefficient becomes negative, is seen in our proof-of-principle case. We show the results of this test in the bottom left panel of Fig. S10, and the DM coefficient indeed prefers a negative value, albeit not to the same degree seen in the real data. In the proof-of-principle example, the DM component prefers a negative flux of magnitude similar to the GCE (rather than several times larger as in the real data); this is consistent with the results of the injection test shown in the main text.

Finally, we can also test the degree to which the DM coefficient floats negative when we simulate data using the same templates employed in the fit; the results are shown in the bottom right panel of Fig. S10. In this case the DM coefficient can still float negative in some realizations, but not to the same extent seen in either the proof-of-principle case or (to a much greater degree) the real data. Specifically, we observe that in 10 out of 20 simulated realizations, the $68\%$ containment bands on the posterior overlap the true DM flux of zero. In the remaining realizations, the upper 68$\%$ containment band on the flux fraction is less than zero, but greater than or equal to $-0.5\%$. In all cases, the degree to which the DM normalization floats negative is representative of how much DM flux needs to be added before DM is recovered in all the analyses.

In all cases, when allowing the DM normalization to float to negative values, all prior ranges remain the same, but the DM normalization prior range is taken to be linear flat, with $A_{\rm NFW}=[-9,9]$.

V Varying diffuse background models

As Galactic diffuse emission is the dominant contribution to the gamma-ray sky and could be a large source of systematic effects, in this section we check that the injection test on the real data behaves similarly to the description in the main text for the alternative diffuse models A and F. These models were identified by Ref. Calore et al. (2015) as providing either the best fit to the data among a wide range of GALPROP-based models (model F) or a comparable fit with physical values for various floated parameters (model A). Both models are expected to provide a better fit to the gamma-ray sky than our default p6v11 diffuse model in the energy band of interest ($2-20$ GeV) Calore et al. (2015).

Table S6 details the prior ranges used for analyses involving diffuse models A and F. A key difference to the previous prior ranges considered in Tab. S1 is that the inverse-Compton and $\pi^{0}$-bremsstrahlung (gas-correlated) contributions are allowed to float separately, whereas the relative normalizations of the $\pi^{0}$ and bremsstrahlung contributions are held fixed in the main analysis with the p6v11 Fermi diffuse model. Note that to speed up the runs, for some parameters we first perform a broad preliminary scan, and then in our final runs we restrict the priors to a narrower range around the converged posterior values. We always use identical priors when directly comparing analyses on real data to those on simulated data.

Tables S7 and S8 show results for the unmasked case for models A and F respectively, both with and without an additional injected DM signal. In the baseline case without injection, the Bayes factor preference for PSs varies strongly between model F, model A, and the p6v11 Galactic diffuse emission model used in the main text, with the smallest value for model F and the largest for the p6v11 model. This suggests that variations in the diffuse modeling can markedly reduce or enhance the apparent significance of a PS population, and consequently that errors in the diffuse modeling could significantly modify the Bayes factors obtained from the real data.

Figure S11 shows the flux posteriors for the unmasked study using diffuse models A and F.

In both Model A and Model F, the injected 2% DM-like signal is excluded by the reconstructed value, as for the p6v11 diffuse model. However, in both cases, the Bayes factor preference for the model with free DM fraction, over the model with the DM component fixed to its injected value, is modest, only $\sim 10$. In this sense model A and model F perform better than the default p6v11 model, in the unmasked case; we note that their Bayes factors in favor of PSs are much smaller than in the p6v11 case, in the data without injection (and much smaller for Model F than Model A).

Tables S9 and S11 show results for the masked case for models A and F respectively. In this scenario model F completely loses any preference for point sources, even without any DM injection. Model A prefers PSs (prior to injection) with a Bayes factor of $\sim 10^{3}$; when extra DM is injected, the Bayes factor in favor of PSs stays fairly constant, but the DM component is not reconstructed correctly.

For diffuse model A we perform simulations for both the masked and unmasked cases, and find that the recovered Bayes factor in the real data is broadly consistent with expectations, albeit on the low side in the masked case (a detailed study of consistency would require a larger number of realizations). For model F, we have performed one simulation in the masked case, and find that there is no strong preference for PSs in the simulated data (consistent with the real data).

Figure S12 shows the flux posteriors for the masked study using diffuse models A and F.

As discussed above, allowing the DM normalization to float negative provides an alternate form of the injection test. The top right panel of Fig. S10 shows the impact of allowing the DM normalization to float negative in model A, when the 3FGL sources are not masked; we find almost identical results for model F. For the unmasked ROI, the DM normalization floats significantly negative in both cases, similarly to the results obtained with the p6v11 Galactic diffuse model.

When the 3FGL sources are masked, and the Galactic diffuse emission is modeled with model A, the DM component also floats negative (consistent with our earlier finding that the injection test is also failed for this combination of model and ROI). When model F is used instead, there is no preference for the DM to float negative, but this is a corollary of the fact that in this case there is no significant preference for the GCE to be comprised of PSs rather than DM.

We note that because we simulate based on the posterior medians, the simulated flux and properties of the NFW PSs are different between different diffuse models; it is thus not straightforward to use these simulations to answer how the sensitivity of the fit to a fixed PS population differs between different diffuse models. It would be interesting to understand whether combinations of diffuse model + ROI that pass the injection test (correctly recovering the injected DM, and with the DM coefficient not floating significantly negative in the absence of injection) retain a robust preference for a dominant PS contribution to the GCE. At present, the example of model F in the masked ROI suggests this is not the case, but this could be either an outlier or simply a reflection of specific features of model F that cause the fit to lose sensitivity to a dominant PS contribution in the masked ROI. We leave an in-depth analysis of this question to future work.

VI Survival function analysis

NPTF methods do not by default reconstruct the likely positions of the PSs. As a proxy, however, we can ask which pixels have photon count values that are particularly unlikely in the absence of PSs – so-called “hot pixels” – following Ref. Lee et al. (2016). Quantitatively, we first perform a Poissonian template fit in the ROI, omitting the PS templates from the standard NPTF pipeline. We then define a reference model as the sum of these Poissonian templates, with normalizations given by the posterior medians from the fit. Suppose the reference model gives an expected mean number of counts $\mu_{p}$ in pixel $p$, and the real data contain $n_{p}$ photon counts in the same pixel. Then we can define the survival function as:

[TABLE]

where CDF$\,[\mu_{p},\,n_{p}]$ is the cumulative probabilty of observing $n_{p}$ counts for a Poisson function with mean $\mu_{p}$. The smaller the value of $\epsilon_{p}$ in a pixel, the more unlikely the observed number of counts are in the absence of PSs, and thus the more likely the pixel is to contain an unresolved PS. Even in the absence of any PSs, there will still be pixels with small $\epsilon_{p}$, purely due to statistical fluctuations; however, the number of pixels with $\epsilon_{p}$ at a given level can discriminate between populations with different photon statistics.

We compute the cumulative fraction of pixels with $\epsilon_{p}<\epsilon$, as a function of $\epsilon$, in the real data; we then repeat this test in simulated data constructed from fits (to the real data) with and without the NFW PS template, and from the Poissonian reference model. Following Ref. Lee et al. (2016), by default we omit disk PSs from both the fits and the resulting simulated data; thus the three cases are (1) Poissonian templates only, (2) Poissonian templates + isotropic PSs, (3) Poissonian templates + isotropic PSs + NFW PSs. For ease of comparison to Ref. Lee et al. (2016), we also use a $1^{\circ}$ radius mask for this analysis, slightly larger than our standard mask; we have tested somewhat smaller masks on a limited number of realizations, and found little effect on the results.

For each case we take the true values of the parameters to be the medians of their posteriors in the relevant fit, and simulate 100 different datasets; the PS populations and Poissonian components are independently redrawn from the source-count functions or template normalizations respectively, for each of the 100 simulations. From each set of 100 simulations, we extract the median fraction of pixels with $\epsilon_{p}<\epsilon$ (as a function of $\epsilon$), and also the standard deviation of this value. We then further repeat this pipeline, but adding a $2\%$ DM flux injection to the real data at all points.

Figure S13 shows the results of this analysis for three Galactic diffuse models: the Fermi p6v11 diffuse model, and two models of Galactic diffuse emission generated using GALPROP Strong and Moskalenko (1998), labeled models A and F in Ref. Calore et al. (2015).

We find that in all cases, the fraction of pixels with small $\epsilon_{p}$ in real data exceeds that for the Poissonian-only reference model, indicating (as expected) the presence of non-Poissonian fluctuations away from the reference model in the data. However, with the p6v11 diffuse model, the simulated data including PSs – either both NFW PSs and isotropic PSs, or isotropic PSs alone – overproduces the fraction of pixels with small $\epsilon_{p}$. This could be due to a mismodeling of the spatial distribution of the PSs, the diffuse background, or both; the pixel $\epsilon_{p}$ depends strongly on where the PS is located, not only its brightness, as even a relatively faint source may be fairly significant if located in a region of low background. For example, suppose a non-isotropic population of PSs lying preferentially in regions of high background brightness is absorbed by the isotropic PS template; these PSs will then be simulated as purely isotropic, without the correlation with bright background regions, and the simulated data will have more PSs lying in regions of low background brightness (and hence more hot pixels) than the real data.

This effect appears to be largely due to the isotropic PS population, which is detected with a significant positive flux in our analysis, unlike the analysis of Ref. Lee et al. (2016) which found it to be consistent with zero. That analysis used an earlier dataset for this CDF test (Pass 7 vs Pass 8), and a slightly different energy range. We have tested changing the energy range in our version of this analysis, to match the results of Ref. Lee et al. (2016), and continue to find that the simulated data overproduces the pixel fraction at small $\epsilon_{p}$ compared to the real data.

We have tested adding disk PSs to the fit and the simulated data, but the fit generally prefers to assign very little flux to disk PSs in our masked ROI, and this change does not qualitatively alter our results. In the other direction, we have tested forcing the isotropic PS population to zero, and fitting purely with NFW PSs; however, we then find the best-fit NFW PS population becomes brighter, with a higher cutoff in flux (presumably because the NFW PS template absorbs sources that were previously attributed to isotropic PSs). In this case the number of small-$\epsilon_{p}$ pixels is again overproduced. However, if we force the NFW PS population to have the same source-count-function parameters as extracted in Ref. Lee et al. (2016) – with the exception of the overall normalization, which is allowed to float – then we recover similar results to Ref. Lee et al. (2016), with the number of small-$\epsilon_{p}$ pixels agreeing fairly well between the real and simulated data.

Replacing the p6v11 diffuse model with models A and F (which give an overall better fit to the data), the mismatch between simulated and real data for small-$\epsilon_{p}$ pixels is significantly mitigated in models A and F. It is interesting to note that in model A, the simulated data containing isotropic PSs but no NFW PSs yields more low-$\epsilon_{p}$ pixels than the case with both isotropic and NFW PSs. This may again be because of the spatial mismodeling issue discussed above; PSs spatially coincident with the GCE, in a high-background region, will generally have larger $\epsilon_{p}$ (lower significance) than PSs with the same brightness distributed isotropically. If the isotropic PS template absorbs such sources, when they are re-simulated, they are likely to be placed in regions with lower background than in the real data, where their significance is accordingly higher.

Using Model A to describe the Galactic diffuse emission, when the artificial DM component is injected into the data, the fraction of pixels with the smallest $\epsilon_{p}$ values remains essentially unchanged. However, the fraction of small-$\epsilon_{p}$ pixels in the simulated data (where the NFW PS template rises to absorb the injected DM component) slightly decreases, becoming more similar to the real data. This is counterintuitive – increasing the number of NFW PSs decreases the fraction of hot pixels – but appears to be due to a shift in the median source-count function for the NFW PSs toward lower fluxes, although this shift is within the uncertainties. This highlights that the error bars for simulated data, in these survival-function analyses, do not include the scatter from the uncertainty in the simulated template parameters.

VII Two Dimensional Posterior Values of Floated Parameters

In this section we show the triangle plots for several key analyses. These are the one and two dimensional posteriors of the parameters floated in the analyses.

Figure S14 shows the triangle plot corresponding to the scenario where we have simulated NP Bubbles, NP Disk, and Poissonian Bubbles, DM, Diffuse and Isotropic backgrounds, when they are analyzed with the same templates but the Bubbles PS are interchanged with NFW PS. The DM flux is not recovered.

Figure S15 shows the triangle plot corresponding to the scenario where we have simulated NP Bubbles, NP Disk, NP isotropic, and Poissonian Bubbles, DM, Diffuse and Isotropic backgrounds, when they are analyzed unmasked with the same templates but the Bubbles PS are interchanged with NFW PS. The DM flux is not recovered.

Figure S16 shows the triangle plot corresponding to the case where the Fermi data are analyzed without any injected DM flux, in the Inner Galaxy ROI. These posteriors for the templates are used as parameter values in the simulated data, except for the DM normalization which is taken from the same run but without NFW PS.

Figure S17 shows the triangle plot corresponding to the analysis of an artificial DM signal inserted into the Fermi data (1.8$\%$ DM injection). The DM flux is not recovered.

Figure S18 shows the triangle plot corresponding to the scenario an artificial DM signal is inserted into the Fermi data (2.4$\%$ DM injection), but with the 3FGL unmasked. The DM flux is not recovered.

Figure S19 shows the triangle plot corresponding to the scenario a larger (6.7$\%$) artificial DM signal inserted into the Fermi data. The DM flux is again not recovered.

Figure S20 shows the triangle plot corresponding to the scenario an even larger (15.2$\%$) artificial DM signal inserted into the Fermi data. Some DM flux is now recovered, but some of it is still misattributed to PS.

We note that in some cases, the index describing the slope of the source count function above the highest break ($n_{1}$) is not well converged within the prior range. We have tried increasing the prior range and do eventually find convergence; however, the physical difference between slopes at the maximum of the prior range and higher values are negligible, both corresponding to an extremely sharp cutoff in the source count function above the highest break.