Gamma-Ray Lines in 15 Years of Fermi-LAT Data: New Constraints on Higgs Portal Dark Matter

Monoenergetic $\gamma$-ray spectral lines are among the cleanest signatures of dark matter annihilation. We analyze 15 years of Fermi-LAT data, find no spectral lines, and place strong constraints on dark matter annihilation to monoenergetic $\gamma$-rays. Additionally, we produce the first double-line analysis of the coupled signals from $\gamma\gamma$ and $Z \gamma$ lines, which proves particularly powerful for dark matter masses above $\sim150$~GeV. From our constraints on a double-line feature, we investigate and constrain some minimal models where the Galactic Center Excess (GCE) can be fit by dark matter annihilation through the Higgs boson into Standard Model particles.


I. INTRODUCTION
Models that generate the observed dark matter (DM) abundance through thermal freezeout provide one of the most compelling explanations for the cosmological evolution of our universe [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Fortunately, scenarios dominated by a 2 → 2 process (e.g. χχ → SM SM) also provide us with a precise, testable target for DM annihilation searches [20,21]. The thermal freeze-out mechanism does not generically predict the Standard Model (SM) final states or branching ratios, and thus it is common to examine DM models dominated by tree-level annihilations to different standard-model particle states, such as bb, τ + τ − , W + W − or other leptonic and hadronic pairs. However, these models are simplified for two reasons. First, DM annihilation may include tree-level couplings to a number of final states, with branching ratios that depend on the decay widths of the intermediate particles. Second, in addition to tree-level processes, there are guaranteed loop-level processes. Some of these final states, like those that produce γγ or Zγ lines, may be more detectable than tree-level annihilation processes despite their subdominant branching fractions.
While the two-photon channel leads to a mono-energetic line at E γ = m DM . The Zγ channel is kinematically accessible at DM masses of m DM > m Z /2, and leads to final state photons with energies centered around The spectrum has an intrinsic width due to the finite life-time of the Z boson, which is given by [22,23] This width leads to an effectively monochromatic signal at the current Fermi-LAT energy resolution. Making use of this fact, we propose a double-line search that further increases our experimental sensitivity. * Email: pedro.delatorreluque@fysik.su.se,ORCID: 0000-0002-4150-2539 † Email: juri.smirnov@liverpool.ac.uk; ORCID: 0000-0002-3082-0929 ‡ Email: linden@fysik.su.se,ORCID: 0000-0001-9888-0971 10 100 10 -29 10 -28

-27
Gamma-ray energy [GeV] 〈σv〉 γγ [cm 3 /s] The relationship between the branching ratios to all final states depends on the mediator choice. The simplest, and most predictive scenario is DM coupling though the Higgs portal, a singlet operator H † H [26][27][28][29][30][31][32][33]. In this case, the branching ratios to all SM final states are entirely fixed by the well-known properties of the Higgs boson. Thus, combining collider-grade accuracy with the freeze-out condition entirely fixes our signal expectation for a given DM mass.
In this work, we reanalyse existing Fermi-LAT data, choosing CLEAN events from 180 months of the PASS8 data [34], and perform both single-and double-line searches for annihilating DM. We find that double-line analyses have a superior constraining power, especially at large DM masses. Furthermore, we show that our limits on Higgs-mediated DM annihilation are in tension with Higgs-portal interpretations of the GCE at masses near the m H /2 resonance [33,35]. Fig. 1 shows the limits of our single-and double-line analyses on the ⟨σv⟩ γγ annihilation cross section, and compares our results with previous work [24,25]. Our analysis provides stronger constraints, particularly at large γ-ray energies.

A. Gamma-ray datasets
The Fermi-LAT is a pair-conversion telescope that measures γ-rays with energies between ∼20 MeV and ∼1 TeV [34]. In this paper, we use ∼180 months of data spanning from 2008-08-04 to 2023-07-20 selecting CLEAN events from the PASS8 data. We include events from good quality time intervals and remove periods when the LAT was operating at rocking angles θ r > 52 deg ((DATA QUAL> 0) && (LAT CONFIG==1) && ABS(ROCK ANGLE)< 52). We also apply the zenith-angle cut θ z < 90 deg, to avoid contamination from the Earth limb. We limit our analysis to EDISP3 events evtype = 512), which have the best energy reconstruction (hence, best energy resolution), to minimize uncertainties relating to instrumental energy dispersion. We employ the P8R3 CLEAN V3 version of the instrument response functions. The extraction of Fermi-LAT data and calculation of exposure maps is performed using the most up-to-date version of the ScienceTools [36] (2.0.8; released on 01/20/2021). In addition, we performed different consistency checks of our analyses using PASS7 and PASS8 front-and back-converted events, which allows us to compare our results with those obtained by the Fermi-LAT Collaboration [37].
We extract data spanning from 10 GeV to 300 GeV in 130 logarithmically uniform energy bins, which constitute an energy resolution ∆E/E ∼ 3% (which is better by at least a factor of two than the intrinsic energy resolution of EDISP3 events). Closely following previous Fermi analyses, we divide our data into three regions of interest (ROIs), which are spherical regions of 3, 16, and 41 degrees around the Galactic center, focusing on the inner regions studied by the Fermi collaboration in Refs. [24,37] (see also [38] for a similar approach). For each ROI, the regions of 68% containment around known sources are subtracted in every dataset following the 4FGL DR2 catalog, except for the ROI3 region, where no source subtraction is done. We also mask the galactic plane in regions with |b| < 5 • and |l| > 6 • as in Refs. [24,37].

B. Single and double line analyses
We closely follow the strategy used by the Fermi collaboration in Ref. [37] and search for spectral lines by performing maximum likelihood fits in each of our ROIs and in 88 sliding energy intervals (the sliding window technique from Ref. [39]) from 10 GeV to 300 GeV. For each interval, we fit the count spectrum in an energy window that surrounds the central energy with a width of ∆E ≥ 3 × σ E (E) 1 (where σ E (E) is the half width of the 68% exposure-weighted energy resolution of each dataset (See [40]), assuming that the background is described by a power-law and adding a linelike signal of free amplitude, which is smeared due to the energy resolution of the Fermi-LAT (Following [41]). We use a likelihood function described by a Poisson distribution in the number of events at every energy window: where N i is the observed number of events at energy E i in each dataset (ROI), and n i is the expected number of counts at energy E ′ that are reconstructed at energy E by the instrument. Under the null hypothesis, there is no line-like signal and the expected number of counts is found by fitting the data to a power-law describing the background emission (n i = n bkg ). In the alternative hypothesis, the number of counts is described as n i = n sig D ef f + n bkg , where D ef f is the energy dispersion matrix that allow us to account for the energy reconstruction of the events by the LAT, and which was obtained using the gtdrm Fermitool.
It is important to remark that, within the energy range where we perform this analysis, the systematic uncertainties in the spectral reconstruction are expected to be negligible compared to the statistical uncertainties in the photon count, as reported by Refs. [24,37]. This assumption does not hold at much lower energies, which would require more complex modeling of Fermi-LAT responses, as in Ref. [24]. The inclusion of systematic uncertainties would only slightly weaken our bounds at low energies, leaving the main conclusions of this manuscript unchanged.
To perform the fits we rely on the Markov Chain Monte Carlo (MCMC) package Emcee [42], since this technique is more robust than conventional optimizers and less prone to finding false local minima. This analysis produces probability distribution functions for every parameter in the fit, which are used to estimate the credible intervals of each parameter and the DM limits. Since the best-fit number of signal events that we obtain is very low, we use the Feldman-Cousins (FC) method [43] 2 to ensure that we are not mis-evaluating the confidence intervals and, hence, the limits. Using the FC method produces roughly the same upper-limits as the MCMC algorithm, except when the best-fit number of source counts is very small. Concretely, we take the best-fit values for the number of background events and number of signal events obtained from the MCMC procedure, and apply the FC method with the likelihood functions defined in Eq. 3. In this way, we reject unrealistically strong upper-limits, particularly for the downward fluctuations.
For the double-line analysis, we repeat the same procedure as the single-line analysis, but add a correlated second line signal that accounts for DM annihilation through the Higgs into a second Zγ line (χ+χ → H → Z +γ). The relative amplitude of the γγ and Zγ signals is correlated by the branching ratio of the Higgs boson to each channel. Moreover, the energy of the line signal for Zγ production is connected to the energy of the γγ signal as described by Eq. 1, and we set E γγ = m χ for the annihilation process that we are considering.
To account for the fact that the energy window in our double line analysis must accommodate both the γγ and Zγ line energies, the lower edge of the energy window considered in this analysis is set to be the minimum value between E Zγ and E γγ − 3σ E , to cover the double-line feature. This results in a larger energy window for the double-line analysis from 50 to ∼ 80 GeV, above which the lower limit of the sliding energy window coincides with the one used in the single-line analysis (i.e. the lower limit is always the 3σ E above 80 GeV).

III. DM BOUNDS FROM THE LINE SEARCH
The expected γ-ray flux from the annihilation of DM particle through the process χ + χ → γ + γ in a region of the sky with angular size ∆Ω is where J ROI (∆Ω) is the astrophysical J-factor that describes the expected annihilation rate given a specific choice of ROI and a DM distribution, dN dE γγ = 2×δ(E −E γγ ) is the γ-ray yield per annihilation, m χ = E γγ is the mass of the WIMP, and ⟨σv⟩ γγ is the annihilation rate to the γγ channel and is related to the total annihilation rate via the mediator-dependent branching fraction ⟨σv⟩ γγ = BR γγ ×⟨σv⟩ ann , where ⟨σv⟩ ann is the total DM annihilation rate.
For the Zγ process, this same formula holds, but with ⟨σv⟩ Zγ = BR Zγ × ⟨σv⟩ ann and dN dE Zγ = δ(E − E Zγ ), where the photon is produced at the energy given by Eq. 1. The main uncertainty in deriving limits on the annihilation rate is the J-factor, J ROI (∆Ω), which directly depends on the Milky Way DM distribution. Here, we assume a local DM density of 0.4 GeV cm −3 [44] and a distance from the Solar System to the GC of 8.5 kpc. We characterize two DM distributions, the NFW and a contracted-NFW profile with an index γ = 1.3 (motivated by studies of the GCE [45,46]), both with a scale radius of r s = 20 kpc.
A. Significance for lines in the gamma-ray spectrum Figure 2 shows the test-statistic (TS) computed for the single-line and double-line analyses as a function of the DM mass. This produces an accurate calculation (assuming Wilks' theorem holds) for the local significance of any line signal (σ local ∼ n line √ n bck ), which can be calculated as Although, the J-factor constitutes the largest uncertainty on the expected annihilation signal from the GC region, we remind the reader that the TS is independent of the J-factor employed. We note no statistically significant peaks (exceeding a 3σ local significance) in any dataset. The most statistically significant peaks hardly exceed 2σ and are not repeatedly present in the different ROIs (i.e. a fluctuation not present in all the datasets). We have performed an analogous analysis, without fixing the branching ratio between the Zγ and γγ processes to the values predicted by the Higgs portal. Also in this general case no significant excess signals were observed and the local significance is roughly identical to the one obtained in the double-line analysis.

B. DM bounds from the single-line and double-line analyses
Given the lack of significant excesses in the γ-ray spectrum, we derive the confidence limits for both the singleand double-line analyses. We produced the bands depicting the 68% and 95% confidence intervals of the observed limits by generating mock data following a power-law distribution with spectral index of −2 and Poissonian noise. We repeat the analysis for 750 iterations (we find that the bands remain stable above ∼ 400 iterations), similar to the approach of Refs. [24,37]. Figure 3 shows the derived limits for ROI41, which provide slightly stronger (but similar) limits than the other ROIs. We show results for other ROIs in the Figure 7 in Appendix C. We also include the observed limits for the double-line analysis as a dashed line, finding that these constraints become stronger at higher energies when the Zγ cross-section becomes larger than the γγ cross section. The ratio BR Zγ /BR γγ is derived from [47] and is discussed in the supplemental material.
In Figure 1, we compare our single-line limits with those obtained by the Fermi-LAT collaboration (using 5.8 yr of data) [24] and those obtained by Foster et al. [48] (using ∼14 years of Fermi-LAT data). The main differences between our analysis and Ref. [48]  a constant energy-window size (of ∆E/E ∼ 0.64), while our analysis utilizes a variable window size based on the local Fermi-LAT energy resolution, (2) their analysis employs a single ROI focused on the inner 30 • around the GC and utilizes a different Galactic plane cut, (3) they utilize SOURCE class photon events with energy reconstructions spanning EDISP 1-3 (the top 75% of well-reconstructed energy events) while we use only CLEAN class photon events from EDISP3 (the 25% of reconstructed events with the best energy resolution) (4) they do not subtract regions surrounding bright γ-ray point sources, while we eliminate these background-dominated regions from our analysis in all ROIs except for ROI3. Despite these differences, our single-line results are in good agreement, as seen in Fig. 1. The DAMPE collaboration recently published the results of their single-line analysis using 5 yr of data collected by the DAMPE instrument [49], obtaining similar bounds to those found in Ref. [24] (see also Ref. [50]).

IV. IMPLICATION FOR HIGGS-MEDIATED ANNIHILATION
Higgs mediated annihilation is unique from the perspective that collider level precision can be used in a DM framework. Since SM processes govern the branching ratios for the annihilation final states, and the DM coupling to the Higgs is fixed by the freeze-out condition, the only unknown parameter is the DM mass.
However, there are several model choices that affect the relationship between the annihilation cross-section and the expected event rates in direct detection and collider experiments [35,51]. Here we mention two well-motivated models: • A singlet scalar model S, with mass m S , which after the electroweak symmetry breaking has the following relevant coupling to the Higgs boson, h where λ p is a dimensionless coupling, and v H is the Higgs field vacuum expectation value [26][27][28]. In this case the spin-independent direct detection cross section is not suppressed, and only λ p < 10 −3 values are compatible with the current limits [32,52]. The GCE signal can thus only be explained in a very narrow mass range around the Higgs resonance [33].
• A Majorana fermion model χ, with mass m χ and the following low-energy couplings where y p is a dimensionless coupling parameter. In this case the direct detection cross section is suppressed for real values of y p [25,35], which means that we would not expect a signal in direct detection searches. Therefore, this model is not constrained by direct detection experiments, and only visible in collider and indirect detection searches.
We note that both scenarios are subject to invisible Higgs decay constraints, which limit BR(h → inv.) < 0.11 [53], and difficult to avoid.
The total annihilation cross section in both scenarios is given by where s is the total s-channel four-momentum, m DM and λ DM are the DM mass and coupling strength to the Higgs, and Γ h is the total Higgs boson decay width. As discussed in Ref. [54] the thermally averaged cross-section at a resonance needs to be performed without the non-relativistic expansion of the annihilation cross section. The fact that the thermal average of the annihilation cross section in the early universe is significantly different from the average at late times leads to a strong late-time mass dependence of the annihilation cross section around the resonance, as discussed in detail in Ref. [32]. Figure 4 shows the Higgs portal parameter space with the predicted total annihilation rate as a function of the DM mass. It is intriguing that γ-ray line searches have the strongest sensitivity around the Higgs resonance, a region in parameter space that is typically challenging to test. We make use of the fact that, given the Fermi-LAT energy resolution, the line signatures can be well distinguished from the continuum emission, such as final state radiation, up to γ-ray energies of E γγ ∼ 300 GeV, as discussed in Refs. [31,33].

V. IMPLICATIONS FOR THE GALACTIC CENTER EXCESS
Observations of γ-ray emission from the Milky Way galactic center have long observed a γ-ray excess that has been The Higgs portal total annihilation rate as a function of the DM mass. We superimpose our constraint from the double-line analysis on the total annihilation rate (red solid). Additionally we show the dwarf spheroidal limits from the Fermi-LAT collaboration [55] (green solid), as well as the constraints from invisible Higgs decay searches [53] (gray dashed). The rate factor predicted by the relic density constraint is shown as a function of DM mass (black dashed).
Since this tree-level final state involves charged particles, there is unavoidably a loop process leading to monoenergetic γγ and, as dictated by electroweak symmetry, narrow Zγ photon lines. However, the simple b−quark loop produces a branching ratio of the order of BR γγ ∼ α 2 EM /(4π) ∼ 10 −6 , which falls far below current experimental sensitivities.
On the other hand, a dominant branching fraction into b−quark states in this mass range is hard to explain unless the interaction is related to the quark Yukawa couplings. Thus, we are naturally led to the Higgs-boson mediated scenario. Notably, a Higgs-motivated bb−annihilation rate that fits the GCE unambiguously predicts bright γγ and Zγ signals. Figure 5 shows the 95% confidence interval for the GCE signal-predictions for the annihilation rates ⟨σv⟩ γγ and ⟨σv⟩ Zγ , as a red and magenta ellipse, respectively. Those regions are derived from the best-fit region of Ref. [60]. We compare the predicted rates with our best limits from the double-line analysis of 15 years of Fermi-LAT data, and find that are results are beginning to be in tension with minimal Higgs portal mediated scenarios as an explanation for the GCE. By minimal we mean models based on the interactions in Eqs. 6, and 7, in which the relic density is determined by the s-channel Higgs annihilation (shown as black dashed line), and given the searches for invisible Higgs decays apply. Note The predictions for the γγ (dark red) and Zγ (magenta) annihilation rates in the Higgs portal model, given the GCE 95% confidence interval found in Ref. [60]. We superimpose our doubleline search limits from a gNFW profile in ROI41 (red solid), and the ⟨σv⟩γγ values predicted by the freeze-out (black dashed). Note that in contrast to constraints from dwarf galaxies this search for γ−lines in the GC does not have an intrinsic J-factor uncertainty, when the comparison to the continuum excess is made. that other approaches find best-fit parameter regions [76], that are only consistent with the relic density predictions within the excluded parameter range.

VI. DISCUSSION AND CONCLUSION
In this letter, we have reanalyzed 15 years of Fermi-LAT data and studied spectral signatures stemming from DM annihliation to monoenergetic lines. We find no evidence for any statistically significant excesses and set strong limits on DM annihilation to mono-energetic photons. Our results can be applied to a broad range of DM scenarios, constraining DM masses up to (or beyond) the electroweak scale in many well-motivated DM models. If non-perturbative effects, such as bound-state formation are taken into account [77][78][79], the reach extends to even higher DM masses.
Additionally, we performed the first double-line analysis of the full Fermi-LAT data-set, using 15 years of data. We placed strong limits on thermal DM that is coupled to the Standard Model through the Higgs portal, in a largely model independent way. Our results are in moderate tension (but do not entirely rule out), Higgs Portal models of the GCE with dark matter masses that sit near the m H /2 resonance. This parameter space is of interest due to the fact that it is the only portion of the Higgs Portal parameter space that is consistent with the GCE and constraints from the branching ratios to invisible particles. Moreover, we note that our constraints are roughly independent of the dark matter density profile near the galactic center, as the cross-sections for both the GCE continuum and the line search shift in the same way.
We emphasise, that the double-line technique can be applied to other datasets, and is particularly promising at energies near the Higgs resonance. Furthermore, it has an enhanced discovery potential for γ-ray signals that have a limited sensitivity due to low photon counts, as it makes use of additional information from photons in correlated energy bins. In this appendix we show examples of the derived count spectra, and the template functions used in the single-and doubleline analyses in the first figure panel. We furthermore show the limits obtained from the single-line analysis for all ROIs studied in this work in the second figure panel, and for the double-line analysis in the third figure panel. Finally, we report the ratio BR Zγ /BR γγ as function of the DM mass, as derived from [47] (consistent with analytic calculations from [80]), which motivates the application of the double-line analysis for O(TeV) gamma-ray data Figure 6 shows the count spectrum for the ROI16 region at two different energies: at around 106 GeV (upper row), where fluctuations at the level of 1 − 2σ are observed for all ROIs (see Fig. 2) and at 155 GeV, where the signals from the photons produced in the γγ and Zγ decay start to merge into the same energy bin (leading to a significantly stronger constraint compared to the one from the single-line analysis above this energy). Here, we include the fitted count spectra assuming only background (null fit) and assuming background+signal (signal fit), which allows us to see a comparison of the signals searched in the singleline (right panels) and double-line (left panels) analysis.
FIG. 6. Count spectra for two different energy windows comparing the double-line (left column) and single-line (right column) best-fit signals.
In the upper panels, the energy window is centered at Eγγ ∼ 106 GeV (EZγ ∼ 87 GeV) and in the lower panels at Eγγ ∼ 155 GeV (EZγ ∼ 142 GeV). At energies greater than Eγγ ∼ 180 GeV both signals expected in the double-line analysis merge.
Appendix B: Single-line limits Figure 7 shows the limits obtained from the single-line analysis compared to those from Refs. [24,48,49]. The ROI16 region constitutes the region that leads to better bounds. Thus, in the upper panels, we show the limits derived assuming an NFW profile, on the left, and a contracted-NFW profile (with index γ = 1.3), on the right. The contracted-NFW profile has been found in several papers to be the DM profile most compatible with the GCE. In this case, the limits from other works have been rescaled accordingly. Then, in the lower panels we show the limits obtained for the ROI3 (left) and ROI41 (right) regions. As we see, the ROI3 region offers the most conservative limits, which is due to the smaller ROI, which produces many fewer counts and larger Poissonian uncertainties, while the ROI41 region leads to similar but slightly higher bounds than in the ROI16 region. Larger ROIs can have large systematic uncertainties associated to the analysis and a lower expected signal-to-noise ratio for a DM signal. 10 [24,48,49]. The upper panels show the limits for the ROI16 region, assuming an NFW and a contracted-NFW (with index γ = 1.3) DM profile, for the upper-left and upper-right panels, respectively. The limits from other works have been rescaled in the case where we show the limits assuming a c-NFW profile. The lower panels show the limits for the ROI3 (left) and ROI41 (right) regions.
Appendix C: Double-line limits Figure 8 shows the DM bounds and confidence bands obtained from the double-line analysis for the ROI16 (left panel) and ROI41 (right panel), compared to the limits derived from the single-line analyses, drawn assuming a c-NFW profile (with index γ = 1.3). This show that the double-line analysis would constitute a more sensitive way to look for line-like gamma-ray signals at very high energies. The reason is that at high energies the photons from the γγ and Zγ decays are produced at roughly the same energy and the branching ratio for the Zγ decay becomes much higher than for the γγ process, leading to a more constraining limit of ⟨σv⟩ γγ . This means that a search for lines at the energies where the Zγ decay dominates lead to a much stronger constraint due to the sum of the signals coming from the γγ and Zγ decays and the branching ratio for the Zγ process.  Figure 9 displays the signal strength ratio between the Zγ and γγ mono-energetic photon signals for Higgs mediated annihilation. At low energies the massive Z-boson emission suppresses the signal, in the intermediate mass range the ratio is determined by the coupling ratio between the weak and electromagnetic couplings α EW ∼ 1/29 and α EM ∼ 1/137, as well as the photon multiplicity. Finally, at larger masses the loop factor of the γγ signal is suppressed due to the destructive interference of virtual particles, which enhances the relative Zγ signal strength.  9. The signal-strength ratio of the Zγ and γγ mono-energetic photon signals for Higgs mediated annihilation [47].