Higgsino Dark Matter Confronts 14 years of Fermi Gamma Ray Data

Thermal higgsino dark matter (DM), with mass around 1 TeV, is a well-motivated, minimal DM scenario that arises in supersymmetric extensions of the Standard Model. Higgsinos may naturally be the lightest superpartners in Split-supersymmetry models that decouple the scalar superpartners while keeping higgsinos and gauginos close to the TeV scale. Higgsino DM may annihilate today to give continuum gamma-ray emission at energies less than a TeV in addition to a line-like signature at energies equal to the mass. Previous searches for higgsino DM, for example with the H.E.S.S. gamma-ray telescope, have not reached the necessary sensitivity to probe the higgsino annihilation cross-section. In this work we make use of 14 years of $\textit{Fermi}$ gamma-ray data at energies above $\sim$10 GeV to search for the continuum emission near the Galactic Center from higgsino annihilation. We interpret our results using DM profiles from Milky Way analogue galaxies in the FIRE-2 hydrodynamic cosmological simulations. We set the strongest constraints to-date on higgsino-like DM. Our results show a mild, $\sim$2$\sigma$ preference for higgsino DM with a mass near the thermal higgsino mass and, depending on the DM density profile, the expected cross-section.

Thermal higgsino dark matter (DM), with mass around 1 TeV, is a well-motivated, minimal DM scenario that arises in supersymmetric extensions of the Standard Model. Higgsinos may naturally be the lightest superpartners in Split-supersymmetry models that decouple the scalar superpartners while keeping higgsinos and gauginos close to the TeV scale. Higgsino DM may annihilate today to give continuum gamma-ray emission at energies less than a TeV in addition to a line-like signature at energies equal to the mass. Previous searches for higgsino DM, for example with the H.E.S.S. gamma-ray telescope, have not reached the necessary sensitivity to probe the higgsino annihilation cross-section. In this work we make use of 14 years of Fermi gamma-ray data at energies above ∼10 GeV to search for the continuum emission near the Galactic Center from higgsino annihilation. We interpret our results using DM profiles from Milky Way analogue galaxies in the FIRE-2 hydrodynamic cosmological simulations. We set the strongest constraints to-date on higgsino-like DM. Our results show a mild, ∼2σ preference for higgsino DM with a mass near the thermal higgsino mass and, depending on the DM density profile, the expected cross-section.
Dark matter (DM) makes up ∼27% of the energy in our Universe today [1], with only ∼5% of the energy density in ordinary matter, yet its microscopic nature remains unknown. One tantalizing possibility, which has driven decades of experimental and theoretical effort, is that DM arises as the lightest superpartner (LSP) in supersymmetric (SUSY) extensions of the Standard Model that address the hierarchy problem related to the unnaturally low Higgs mass parameter (see [2] for a review). LSP DM at the ∼TeV scale may naturally acquire the correct DM abundance through thermal freeze-out. On the other hand, electroweak scale SUSY and LSP DM have come under increasing tension in recent years from null searches for new physics at the Large Hadron Collider (LHC) [3,4], direct detection experiments [5], and indirect searches [6,7]. Natural LSP candidates are the neutral gauginosnamely, bino and wino LSPs -and the higgsino. Pure bino DM is excluded by direct searches at the LHC [8,9]. Wino DM is in strong tension with null results from DM annihilation searches with the H.E.S.S. gamma-ray telescope [6,10]. Nearly-pure higgsino LSPs, on the other hand, remain one of the better motivated and sought after, yet unprobed, DM scenarios (see, e.g., [11] for a recent summary). Higgsinos, which are the superpartners of the two Higgs doublets in the minimal supersymmetric standard model (MSSM), are especially motivated in light of (i) null results for wino DM, and (ii) the fact that null searches for superpartners at e.g. the LHC suggest that nature may implement a split-spectrum version of SUSY such as Split-SUSY [12][13][14][15][16], which naturally leads to higgsino or wino LSP DM. Split-SUSY, mini-Split, and similar constructions [17] aim to preserve LSP DM and high-scale gauge unification but give up on trying to fully solve the hierarchy problem; in such models the scalar superpartners are taken to have large masses, with the gauginos and higgsinos remaining near the TeV scale. Such split-spectrum models may accommodate the observed Higgs mass and solve a number of troublesome problems with the MSSM, such as the lack of flavor changing neutral currents [18] and the nonobservation of new CP violation in electric dipole moment searches [19,20].
In this work we search for annihilation signatures of higgsino DM with Fermi gamma-ray data. The higgsino interactions with the Standard Model are specified by its representation under the electroweak force, and thus the requirement that freeze-out of annihilation produces the correct DM abundance fully determines the higgsino mass under thermal cosmologies. The higgsino mass arises from the Lagrangian term L ⊃ −µH u · ·H d + h.c., where µ is the MSSM µparameter,H u (H d ) is the up-type (down-type) higgsino electroweak doublet, and is the totally anti-symmetric symbol in SU (2) L space. There are two neutral higgsino fermions, which are generically split into two nondegenerate Majorana mass eigenstates by dimension-five operators that have the effect of inducing a slight mixing between the neutral gauginos and higgsinos (see, e.g., [21]). The charged higgsino states are heavier than the neutral states by at least ∼350 MeV because of radiative contributions to the charged higgsino masses below electroweak symmetry breaking [22]. The mass splitting between the two neutral Majorana states, which we call ∆m, must be greater than around 200 keV to avoid direct detection constraints from inelastic Z-exchange, where the lower Majorana state scatters into the heavier mass eigenstate (see, e.g., [21]). When ∆m 200 keV, direct detection of higgsino DM proceeds through elastic scattering with higher-dimensional operators and is thought to be below the neutrino floor [23][24][25][26].
The relic abundance of higgsinos from thermal freezeout matches the observed DM abundance [1] for mass m χ = 1.08±0.02 TeV, accounting for uncertainties on the DM abundance [27]. We refer to the higgsino with such a mass as the thermal higgsino. Apart from the SUSY motivations, higgsino DM may be viewed through the lens of minimal DM [28]. If m χ is less than the TeV scale, then higgsinos are a computable but subdominant component of the DM, unless there is a non-standard cosmological history that increases their abundance. The same annihilation processes that set the higgsino DM abundance in the early universe also lead to annihilation signatures of higgsinos today. Sommerfeld enhancement introduces dependencies of the present-day annihilation rate on the mass splittings between the neutral and charged higgsino states, though this effect is relatively minor at the thermal higgsino mass [29].
The strongest existing indirect detection constraints on the thermal higgsino arise from Galactic Center (GC) searches for the line emission expected due to the γγ and γZ final states with H.E.S.S. [30], which constrain σv γγ 4 × 10 −28 cm −3 /s assuming an Einasto DM profile, whereas the thermal cross section is a factor of 4 smaller, and for H.E.S.S. searches for continuum emission from annihilation to W + W − [31]. The forthcoming Cherenkov Telescope Array (CTA), on the other hand, with 500 hours of exposure, is expected to have sensitivity to higgsino DM at the thermal mass [32]. Future lepton or hadron colliders may also be able to discover thermal higgsinos [33][34][35].
In this Letter we use the existing 14 years of Fermi data to achieve world-leading sensitivity to higgsino DM by searching for annihilation to continuum gamma-rays through the W + W − and ZZ final states. Our upper limits on the annihilation cross-section surpass those from H.E.S.S. for higgsino-like DM with mass m χ ∼ TeV, though our upper limits are weaker than expected due to the presence of a modest (∼2σ) preference for the signal model over the null hypothesis. We interpret our results in the context of DM profiles from the FIRE-2 hydrodynamic cosmological simulations [36, 37] to show that the best-fit annihilation cross-section we recover may be consistent with the expected higgsino cross-section, potentially providing the first hint of thermal higgsino DM. Data reduction and analysis.-We reduce 722 weeks of Pass 8 Fermi gamma-ray data with SOURCE selection criterion taken between August 4, 2008 and June 10, 2022 with the recommended quality cuts DATAQUAL>0 and LATCONFIG==1 along with zenithangle less than 90 • . We include the top 3 of 4 quartiles of the data as ranked by the point spread function (PSF). As in [38], we initially bin the data into 40 logarithmically-spaced energy bins between 200 MeV and 2 TeV, and we bin spatially using HEALPIX [39] with nside=512. However, in our analysis we only analyze data starting in the 18 th energy bin, with minimum energy ∼ 10.02 GeV, since our signal peaks at higher energies and since the lower energies are  Figure 1. The photon count data used in this work in our ROI, which has the Galactic plane masked at |b| ≥ 1 • along with 4FGL PSs and pixels more than 10 • from the GC. For illustration the data are summed above 10 GeV. We analyze the data in 9 concentric annuli, as indicated.  more contaminated by Galactic diffuse emission. Starting at 10 GeV we also mostly avoid the Fermi Galactic Center Excess, which is an excess of ∼GeV gamma-rays observed near the GC [40-45]; the excess has not been found to extend above 10 GeV with our Galactic emission model [46]. We include energies up to the DM mass m χ , as the signal spectrum has no support beyond that.
Our region of interest (ROI) for the analysis is that within 10 • of the GC, with the Galactic plane masked (|b| ≥ 1 • ) in addition to a 4FGL point source (PS) mask [47]. The PS mask is constructed through the following procedure. First, we use [47] to construct a PS model in our first analysis energy bin in terms of predicted photon counts accounting for the detector response and, in particular, the PSF. Then, we find the brightest pixel, as ranked by predicted PS emission, within 30 • of the GC. We mask all pixels that have a predicted photon flux from PSs as low as 0.5% of that in the brightest pixel. This masking approach preferentially masks pixels from bright sources, as opposed to more conventional masks based off of PSF containment that may under-(over-) mask bright (dim) sources. We then further divide our ROI into 9 concentric annuli, going out to 10 • from the GC starting at 1 • , with angular spacings of 1 • . We stack and analyze the spectral data in each of these annuli independently. The photon counts in our analysis ROI above 10 GeV are illustrated in the top panel of Fig. 1. Even after the Galactic plane mask, more photons are observed near the GC and the Galactic plane from diffuse emission within the Milky Way.
We model the spectral data in each annulus under the null hypothesis using a linear combination of: (i) the spectral template derived from the Fermi Galactic emission model gll iem v07 (p8r3), reprocessed for our data set and selection criterion; (ii) the 4FGL PS spectral template appropriate for our ROI; and (iii) the isotropic diffuse emission appropriate for our Galactic emission model and data set. PS and isotropic emission, however, are sub-dominant compared to Galactic diffuse emission, as illustrated in Supplementary Material (SM) Fig. S1. As shown in the SM, not including the PS and isotropic templates leads to nearly identical results. We define the ensemble of null-hypothesis spectral models (Galactic emission, PS, and isotropic) as the background templates. In each radial bin we construct a likelihood to constrain the spectral model, which consists of the background templates along with the signal template that is discussed shortly, by taking the product of the Poisson probabilities to observe the data counts in each energy bin given the model prediction. The background templates are given individual nuisance parameters that rescale the overall normalization of that template; we require the nuisance parameters to be positive. At a given DM mass m χ and in a given radial annulus, we construct the profile likelihood for the annihilation cross-section σv profiling over the background nuisance parameters. We then construct the joint profile likelihood for σv , at fixed m χ , by taking the product of the profile likelihoods over all radial bins. We use the joint profile likelihood to constrain the signal model. (See the SM for details, along with an alternative analysis that incorporates spatial information into the likelihood.) In Fig. 2 we show the background-subtracted counts data, with the best-fit null hypothesis model, summed over all annuli up through 1.1 TeV for an analysis looking for a higgsino with m χ = 1.1 TeV. The data are largely consistent with the null hypothesis. Note that this figure is for illustrative purposes only and is not used in the analysis, which treats the radial bins separately. Our sensitivity is dominated by the energy bins less than around 100 GeV, as we further illustrate in the SM; this should be contrasted with the sensitivity of upcoming experiments like CTA that will probe thermal higgsino DM at energies near a TeV but lose sensitivity below ∼100 GeV. Our inclusion of photons with energies between ∼10-100 GeV is what makes us competitive with CTA, even though CTA will have a much larger effective area than Fermi.
Results.-To interpret the data in the context of the higgsino model we need to compute the gamma-ray spectrum dN/dE per annihilation from the decays of the unstable particles produced during higgsino annihilation. The dominant annihilation channels are χχ → W + W − and χχ → ZZ. At m χ = 1.1 TeV the branching ratio to W (Z) pairs is ∼60% (∼40%). Since the higgsino annihilates through its electroweak interactions, for a given m χ there is a fixed and calculable σv . This annihilation cross-section is illustrated in Fig. 3 as a function of m χ . Note that there is minor dependence on the mass splittings between higgsino states, though as we show in the SM these differences do not qualitatively affect our results. The dN/dE for annihilations to W and Z pairs are calculated using PPPC 4 DM ID [48]. When constraining the higgsino model we treat the overall cross-section σv as a free parameter, which could be negative, but with the branching ratio to W and Z pairs fixed.
In addition to the spectrum per annihilation we also need to know the astrophysical J-factor, J ≡ dsρ 2 DM (s), in our ROI in order to compute the expected signal in our radial annuli. Here, ρ DM (s) is the DM density along the line-of-sight, parameterized by the distance s from Earth. Our benchmark DM profile is the spherically-symmetric Navarro-Frenk-White (NFW) [49, 50] profile, normalized to produce a local DM density ρ DM (s = 0) = 0.4 GeV/cm 3 , with a scale radius r s = 20 kpc, and with the distance from the GC to the Sun of r = 8.23 kpc [51]. Our local DM density choice is motivated by the recent review [52], which concludes that local DM density measurements, from analyses of stellar motions perpendicular to the disk within a few kpc, tend to favor the range 0.4 − 0.6 GeV/cm 3 , broadly consistent with rotation curve data, though one should keep in mind that the true local DM density may be slightly larger or smaller than our choice. Furthermore, the scale radius is currently poorly constrained, such that r s = 20 kpc represents a reasonable choice as opposed to a value strongly preferred by data.
Given the ad-hoc nature of the NFW profile, and that it is motivated by DM-only N -body simulations, a potentially more promising approach to computing the Jfactor profiles is to use the results for Milky Way analogue galaxies in hydrodynamic cosmological simulations that include baryonic effects. Towards that end, we also compute the J-factor profiles in 12 FIRE-2 zoom-in Milky Way analogue galaxies [36], using simulation outputs provided in [37]. The FIRE-2 simulations are state-of-theart hydrodynamic simulations, which provide the highest angular resolution to-date for J-factor profiles, that  . The best-fit annihilation cross-section for higgsino-like DM as a function of the DM mass for our fiducial analysis assuming an NFW DM profile (left) and the FIRE-2 Romulus profile (right), which gives the best-fit to the data of all profiles considered. We illustrate the best-fit, 1σ (green) and 2σ (gold) confidence intervals for the recovered cross-section, in addition to the 95% one-sided upper limit. We compare our results to the expected higgsino annihilation cross-section (red) and to the 95% upper limits from the H.E.S.S. searches for annihilation to W + W − and gamma-ray lines. The mχ range, accounting for uncertainties, where higgsinos make up the correct DM abundance in the standard thermal cosmology is shaded (thermal mχ). For the Romulus profile the recovered cross-section is consistent with the expected cross-section for the thermal higgsino within 1σ (green band) and inconsistent with the null hypothesis of no higgsino at ∼2σ, as illustrated by the gold band. Note that σv is allowed to be both positive and negative in the analysis, even though negative cross-sections are unphysical. For all mχ we assume that higgsinos make up all of the observed DM.
incorporate e.g. stellar feedback and radiative transfer amongst baryons, which dominate the inner potential wells of Milky Way sized galaxies, in addition to gravitational dynamics. Six of these twelve galaxies, including Romulus and Romeo which we discuss more below, were evolved in pair configurations to mimic the interactions between the Milky Way and Andromeda. The Milky Way analogues are chosen to have stellar masses in the range (3, 11) × 10 10 M with virial masses in (0.9, 1.8) × 10 12 M [37]. The particle masses and positions were then adjusted in [37] such that the local DM density is 0.38 GeV/cm 3 at the distance to the Sun r = 8.3 kpc, which are similar to our fiducial values for the NFW profile. We then compute the azimuthally-averaged J-factors in our ROI annuli; these J-factors are compared to those from the NFW profile in the SM Fig. S5.
In the left panel of Fig. 3 we show the results of our analysis of the Fermi data interpreted for higgsino DM using the NFW DM profile. We illustrate the best-fit cross-section, along with 1 and 2σ significance containment intervals, as functions of the higgsino mass, assuming that at each mass the higgsino makes up all of the DM (see SM Fig. S10 for our results assuming a subfraction of the DM). Our one-sided 95% upper limit is also illustrated. For the fiducial NFW profile we are unable to exclude the higgsino cross-sections over the mass range shown. However, our upper limit is weaker than expected due to a slight statistical preference for the sig-nal model over the null model. At the thermal mass the local significance in favor of the signal model is ∼2σ (see SM Fig. S8 for the discovery test statistic as a function of mass). Note that in Fig. 2 we illustrate the higgsino model prediction relative to the background-subtracted and fully-stacked data for a reference cross-section.
The FIRE-2 J-factor profiles are typically enhanced relative to that of the NFW model due to adiabatic contraction, as illustrated in SM Fig. S5, though there is significant spread over the 12 realizations. In the right panel of Fig. 3 we show our results interpreted in the context of the FIRE-2 halo, Romulus, with the largest J-factor (the other halos are illustrated in the SM). Note that we use the FIRE-2 naming conventions for the Milky Way analogue galaxies [36]. The FIRE-2 halo profiles lead to comparable discovery significances compared to the NFW analysis, as illustrated in SM Fig. S8, with Romulus providing the best fit. Intriguingly, with the Romulus profile and multiple other FIRE-2 profiles the excess in favor of the signal model has a best-fit crosssection consistent with the higgsino model at the thermal mass. Over the ensemble of 12 FIRE-2 J-factor profiles that we consider, the best-fit σv for m χ ≈ 1.08 TeV ranges from 1.7 × 10 −26 cm 3 /s to 7.3 × 10 −26 cm 3 /s, with the median value of 3.4 × 10 −26 cm 3 /s; the higgsino cross-section at this mass is σv ≈ 1.3 × 10 −26 cm 3 /s. The Romulus profile leads to the best-fit cross-section at m χ ≈ 1.08 TeV of σv = 1.7 ± 0.8 × 10 −26 cm 3 /s. The Romeo halo may be the most Milky Way-like, due to the similarities of its thick disk, circular velocity, and stellar mass to the Milky Way; using this halo we recover a cross-section σv = 1.9 ± 1.0 × 10 −26 cm 3 /s at the thermal higgsino mass.
In Fig. 3 (left) we show the 95% upper limit from an analysis of H.E.S.S. data looking for continuum emission above ∼200 GeV associated with χχ → W + W − [31]. H.E.S.S. is less sensitive to χχ → ZZ, since annihilation to Z pairs produces significantly fewer photons above ∼200 GeV than annihilation to W pairs, so to convert the results presented in [31] to higgsino-like DM limits we use only the W + W − result (additionally, H.E.S.S. does not present results for annihilation to Z pairs). Furthermore, [31] uses an ROI ranging from 0.5 • to 2.9 • from the GC and assumes an Einasto profile; we rescale their results to those appropriate for an NFW profile in the left panel of Fig Constraints on higgsino DM using H.E.S.S. searches for gamma-ray lines, which are from the loop-suppressed processes χχ → γγ and χχ → Zγ, are also relevant, though the hard-photon spectrum is affected by electroweak radiative effects [53]. We translate the H.E.S.S. gamma-ray line limits in [30], which were computed using the same ROI as in their continuum W + W − search described above, to limits on the total annihilation crosssection using the NLL' calculation for the energy spectrum near the gamma-ray endpoint in [53]; the recasted limit is illustrated in the left panel of Fig. 3. The H.E.S.S. upper limit surpasses our upper limit at large masses, though it should be kept in mind that (i) the H.E.S.S. analysis is significantly closer to the GC than ours and thus the comparison relies on the possibly incorrect shape of the NFW profile, and (ii) the result in [53] may be subject to O(1) uncertainties when applied to the H.E.S.S. analysis because the energy binning used in the analysis does not directly match the assumptions in the calculations in [53]. (See also [54].) Discussion.-In this work we set the strongest constraints to-date on higgsino-like DM that annihilates to W + W − and ZZ using nearly the entire Fermi data set collected since the mission's launch in 2008. We search for the continuum gamma-ray emission above 10 GeV associated with the decays of these massive vector bosons. Interpreting our results in the context of the NFW profile our upper limits on the annihilation cross-section are world-leading but do not constrain the predicted higgsino cross-section. On the other hand, hydrodynamic simulations show that the DM profile may be contracted in the inner Galaxy, enhancing the annihilation rate.
Our upper limits are weaker than expected because of a modest (∼2σ) preference for the higgsino model over the null hypothesis of background-only emission. The best-fit cross-section is consistent with the expected higgsino cross-section for a thermal (m χ ≈ 1.1 TeV) higgsino making up all of the DM for multiple FIRE-2 DM density profiles. Given that higgsino DM is well motivated from supersymmetry and, in particular, Split-SUSY type models that seem to be favored by the lack of observed superpartners to-date at the LHC in addition to the observed Higgs mass and precision Grand Unification, the possibility that the data present the first hint for higgsino DM is promising. This possibility will be tested with the upcoming CTA [32], which should be sensitive to higgsino DM annihilation.   Sec. I provides additional details of the methods used in our analyses, Sec. II provides additional results for our fiducial analysis presented in the main Letter, while Sec. III gives the results of systematic analysis variations that provide additional support to our main conclusions.

I. METHODS
In this section we describe the methods used in our fiducial analysis and in our analysis variations.

A. Data analysis
As described in the main Letter, the Fermi data is initially reduced into 40 log-spaced energy bins between 200 MeV and 2 TeV and into spatial maps using the HEALPIX pixelation [39] scheme with nside=512. While reducing the data we also reprocess the p8r3 Galactic emission model, along with the 4FGL PS model, extended source emission model, and the associated isotropic flux model. Note that the p8r3 model is not recommended for use in searches for extended emission because detected, extended emission regions (specifically, unassociated sources with extension beyond 2 • ) were added in as residuals to the p8r3 model. Our signal, corresponding spatially to the J-factor profile, is extended beyond the instrument PSF. However, our fiducial analysis only uses spectral information in the diffuse model and not spatial information. Note also that there are no significant extended sources in the extended source catalog within our ROI, and in addition the GC Excess below 10 GeV has been successfully observed with the p8r3 model [46]. We also use the p8r3 Galactic emission model in our spatial template analysis described below, though ideally it should be verified if Fermi added any residuals into p8r3 that are spatially coincident with our signal template. In future work it would also be useful to perform such spatial template analyses with Galactic emission models that do not include data-driven residuals. Note that later in the SM, in Sec. III D, we also present spectral analysis results using a purely data-driven background model that is constructed from data further away from the GC than our signal ROI.
We perform two types of analyses in this work. First, our fiducial analysis, which is described in the main Letter, uses spectral information only within 9 annuli centered at the GC (see Fig. 1 for the geometry). The second analysis is used as a systematic check in the SM, and it consists of spatial template fits performed independently in each energy bin. We describe both analyses below, starting with our fiducial (spectral) analysis. Note that in all of our analyses we focus on the GC because it provides the most sensitive probe of DM annihilation; searches using dwarf galaxies, for example, would likely be less sensitive [55].
Consider the data d = {n i } stacked over pixels within a given annulus, numbered one through nine with one being the inner-most annulus from 1 • to 2 • . Here, n i stands for the number of counts observed in energy bin i, with i in principle going from 1 to 40 though in practice we restrict the analysis energy range to energies above 10 GeV (for our fiducial analysis) as well as energies below the DM mass. In our fiducial analysis we describe the data using a spectral model that predicts the expected number of counts µ i (θ) in energy bin i as a function of the model parameters θ. In our fiducial analysis θ = {A sig , A bkg }, with A sig the signal model parameter that controls the strength of the signal (i.e., the cross-section), at a fixed DM mass, while A bkg = {A gal , A PS , A iso } changes the normalization of the associated background model components. In particular, µ i (θ) = A sig µ sig i +A gal µ gal i +A PS µ PS i +A iso µ iso i , with µ gal the Galactic diffuse spectrum, µ PS the residual PS spectrum (which extends beyond our PS mask), and µ iso the isotropic spectrum as predicted by Fermi for our event class and diffuse model. The background parameters are not allowed to take on negative values, though our signal parameter is, and the background parameters are also constrained to be less than five times their expected normalizations to aid with convergence of the global optimizer used to maximize the likelihood.
In Fig. S1 we show the spectral templates summed over our first two annuli (1 • < r < 3 • ) with default normalizations (A bkg = 1) and a signal normalization corresponding to σv = 5×10 −26 cm 3 /s for an NFW DM profile and m χ = 1.1 TeV. The Galactic diffuse template dominates over PS and isotropic emission. Indeed, as we show later removing those latter two templates completely does not qualitatively affect our results. Note that we also illustrate an Inverse Compton (IC) spectral template in Fig. S1; this component is not used in our fiducial analysis, but it is used in a systematic test described in Sec. III F.   Figure S1. The spectrum within our first two annuli compared to the spectral templates used in our fiducial analysis, including the p8r3 Galactic emission model, the residual PS model, and the isotropic model. The signal model is illustrated for a DM mass of mχ = 1.1 TeV and a cross-section of σv = 5 × 10 −26 cm 3 /s. Note that we do not include line emission in the signal spectrum, with the peak near a TeV arising from electroweak radiative effects in the W + W − spectrum. We also show the IC spectral template used in Sec. III F. We constrain the model parameters using the data through a Poisson likelihood with the product being over energy bins. In particular, we construct the test statistic (TS) for upper limits where in the second termθ refers to the model parameters that maximize the likelihood while in the first termÂ bkg refers to the background nuisance parameter vector that maximizes the log likelihood at fixed A sig . We assume Wilks' theorem such that the 95% one-sided upper limit is given by the value of A sig greater than the best-fit valueÂ sig where t(A sig ) ≈ 2.71 (see, e.g., [56]). The discovery TS for a one-sided test, which determines the significance of the signal model over the null model, is given by q = t(A sig = 0) ifÂ sig > 0 and q = 0 otherwise. Note that we allow A sig < 0, even though such values are unphysical, to ensure that we find the likelihood maximum. In the SM we perform an alternate analysis using a spatial template fit instead of a spectral template fit. In each energy bin we assign independent nuisance parameters for the PS template and the Galactic emission template. The spatial PS template is constructed from the 4FGL catalog along with the instrument PSF, though we treat the overall normalization of the template as a nuisance parameter. In each energy bin we assign independent nuisance parameters to the PS model and to the Galactic emission model. Thus, the likelihood may be written as where n ij is the observed number of counts in energy bin i and spatial pixel j, while µ ij is the model prediction in that pixel. There is still a single signal parameter A sig , which -given a DM mass -rescales the signal amplitude in each energy bin according to the DM spectrum. The spatial profile of the signal is given by the J-factor. The number of nuisance parameters is twice the number of energy bins.

B. Higgsino annihilation signal
For a nearly-pure higgsino, the freeze-out abundance under a thermal cosmology is essentially determined by its mass, and thus the requirement that it saturates the total amount of observed DM fixes a precise mass prediction.
In the present day, when the higgsino is much colder and Sommerfeld effects become important, the cross-section of its annihilation and branching ratios of its decay products acquire a weak dependence on the two mass splitting parameters: that between its two neutral components and that between the LSP and charged component. In this Letter and associated SM, we present results for the scenario where the charged splitting is fixed to its radiative value 350 MeV and the neutral splitting saturates the up-scattering bound from nuclear recoil experiments, 200 keV [32, 57].
The differential flux of γ-rays expected from higgsino annihilation is given as The astrophysical J-factor encodes information on the DM distribution within the ROI. Note that the higgsino need not constitute the entirety of the observed DM abundance, and indeed the thermal abundance of a higgsino LSP with m χ < 1.1 TeV is a fixed fraction of the observed total, Ω χ /Ω cdm ∼ (m χ /1.1TeV) 2 . However, this fractional abundance is relatively easily enhanced with non-thermal cosmologies.
In this paper we search for only the continuum part of the higgsino annihilation spectra; while this annihilation does produce a monochromatic line signal that is generally easier to target in searches, the expected strength of this signal is comparatively suppressed even after accounting for Sommerfeld enhancement. The continuum spectra consists of contributions from W + W − and ZZ final states (and a small contribution from the γZ channel weighed by a factor of 1/2), and the cross sections are computed at tree-level. These cross sections are modified by the mixing of the various higgsino states due to Sommerfeld enhancement, which is parametrized by where the mixed-state vector v i is determined by solving the Schrödinger equation, whose potential is set by the mass gaps between higgsino states and the electroweak interaction α W e −M r /r [29, 32, 58]. The injection spectra of these annihilation modes into gamma-rays, dN γ /dE, are obtained via the prescription of PPPC 4 DM ID [48], where radiative electroweak corrections have been accounted for and final W/Z polarizations averaged over. Our main analysis assumes a signal that consists of only the prompt emission in order to minimize dependency on propagation functions, and we do not expect secondary contributions to have a significant impact on our analysis. The expected primary flux from a 1.1 TeV thermal higgsino for a fiducial halo profile is shown in Figure S2.

II. EXTENDED RESULTS FOR THE FIDUCIAL ANALYSIS
In this section we present additional results from the fiducial analysis described in the main Letter. In Figs. S3 and S4 we show the best-fit null spectra and residuals in each of the 9 concentric annuli. We illustrate the higgsino signal on top of the background-subtracted data for an NFW DM profile with a higgsino mass m χ = 1.1 TeV and a cross-section σv = 5 × 10 −26 cm 3 /s. Note that these figures are for presentation purposes only and are not used in our analysis, as our analysis profiles over the background nuisance parameter when constructing the profile likelihood for the cross section. However, visual inspection of these figures shows that there are no obvious sources of systematic mismodeling present in our fiducial analysis.
In the main Letter we presented results for the NFW DM profile and for the Romulus FIRE-2 profile, which provided the best-fit to the Fermi data. Here, we present results for all 12 of the FIRE-2 J-factor profiles for Milky Way analogue galaxies. The J-factors for the FIRE-2 Milky Way analogue galaxies are illustrated in Fig. S5 (reproduced from [37]), compared to the J-factor profile for our fiducial NFW profile. Interestingly, almost all of the FIRE-2 J-factors are larger in the inner few degrees relative to those from the NFW profile. On the other hand, the angular resolution in the FIRE-2 simulations, which have the best resolution of any current hydrodynamic simulation suites, is estimated at 2.75 • [37], meaning that the J-factors in the inner ∼3 degrees may be unreliable.
In Figs. S6 and S7 we show the results of our analysis, as in Fig. 3, interpreted in the context of higgsino DM for all 12 FIRE-2 galaxies. Note that we also compare the 95% upper limits to that found using the NFW DM profile. In 11 of 12 FIRE-2 galaxies we find a stronger upper limit than we do with the NFW profile, suggesting that adiabatic contraction with baryons likely enhances the DM signal in Milky Way type galaxies relative to the NFW expectation.
In Fig. S8 we show the discovery TS as a function of the higgsino mass m χ for analyses using the NFW DM profile and the 12 FIRE-2 J-factor profiles. The evidence in favor of the signal hypothesis is ∼1-2σ in significance at all masses. Thus, we conclude that the data shows no evidence for higgsino DM. The TS is slightly larger than expected under the null hypothesis, which could indicate the first hint of a DM signal, but this modest TS could also arise from one or a combination of mismodeling and random statistical chance.
In Fig. S9 we show the best-fit cross-section and 1σ confidence interval in each radial bin independently, interpreted using the NFW DM profile (left panel) and the Romulus profile (right panel), with the higgsino mass fixed to be m χ = 1.1 TeV. These figures show that the analyses in the independent annuli, which are joined together in our fiducial analysis, produce consistent results up to statistical fluctuations.
Throughout the main Letter we assume that regardless of the mass m χ , the higgsino makes up 100% of the DM. An alternative possibility, however, is that for m χ less than the thermal mass the higgsino is a subfraction of the total DM, with the fraction scaling like (m χ /1.08 TeV) 2 for a thermal mass of 1.08 TeV and the standard thermal cosmology. In this case, the DM annihilation signal is reduced by a factor (m χ /1.08 TeV) 4 for m χ < 1.08 TeV. In Fig. S10 we illustrate our results assuming that the higgsino is a subfraction of the DM for a thermal mass of 1.08 TeV. Note that we do not show masses above 1.08 TeV since such masses would overproduce the DM in this framework.   Figure S3. (top panels) The observed flux in the first four annuli, as indicated, along with the best-fit null model that is the p8r3 Galactic emission model plus isotropic and PS emission. (bottom panels) The residuals -data minus best-fit background -constructed from the information in the top panels. We illustrate an example signal for a higgsino with mχ = 1.1 TeV and σv = 5 × 10 26 cm 3 /s for the NFW DM profile.   Figure S5. The J-factor profiles averaged over our analysis annuli for our fiducial NFW profile and for the 12 FIRE-2 J-factor profiles from [37] for Milky Way analogue galaxies. Note that the resolution of the FIRE-2 zoom-in simulations is indicated in shaded grey, which implies that the first 2 annuli may misestimate the true J-factors in those galaxies. The Romeo galaxy may be the most Milky Way like, while we illustrate our results in the main Letter for the NFW and Romulus profiles. The named galaxy simulations, including Remus and Romulus, were evolved in pairs in order to mimic the Milky Way's co-evolution with M31, and are in this sense more realistic representations compared to the six unnamed ones. . We compare the 95% upper limits to that found using the NFW profile; in 11 of 12 cases the FIRE-2 galaxies lead to stronger upper limits relative to that from the NFW profile.   Figure S8. The discovery TS in favor of the signal model over the null hypothesis. The Z-score for discovery is given approximately by the square root of the TS for a one-sided test. We illustrate the TS as a function of mχ for analyses using the NFW DM profile and the 12 FIRE-2 J-factor profiles.   Fig. 3 but instead of assuming the higgsino is all of the DM regardless of mχ we assume that the higgsino is a subfraction of the DM as computed under the assumption of a standard thermal cosmology and with a thermal mass of mχ = 1.08 TeV. We do not show our results for masses above 1.08 TeV in this case since the DM abundance is overproduced in this region given our assumptions for this figure.

III. SYSTEMATIC ANALYSIS VARIATIONS
In this section we describe the results of systematic variations to our fiducial analysis. Before describing the variations in detail, we refer to Fig. S11, which summarizes the expected changes to sensitivity for a subset of the systematic analysis variations described below in Sec. III A and Sec. III B. In Sec. III A we consider changing our analysis ROI, while in Sec. III B we change the analysis energy range. In Fig. S11 we show the expected 95% upper limit in our fiducial analysis assuming the NFW DM profile (left panel) and Romulus profile (right panel). These expected upper limits are computed using the Asimov procedure [56], whereby we assume the data is given by the best-fit null hypothesis model. We also illustrate our fiducial best-fit cross-sections from the analyses on the actual data.
We compare our fiducial analysis expectations to those from the Asimov procedure as applied to our analysis variations in Sec. III A and Sec. III B. In particular, we show the effects of reducing the ROI and changing the minimum analysis energy, as indicated. We illustrate these changes through the expected 95% upper limit, as the changes to the other quantities of interest may be derived from this information. All of our analysis variations have comparable sensitivity, making the cross-checks described in Sec. III A and Sec. III B powerful probes of possible mismodeling.  Figure S11. The best fit annihilation cross-sections for our fiducial analysis compared to the expected 95% upper limit for our fiducial analysis and some of the analysis variations considered in Secs. III A and III B.

A. Reduced ROI
In Fig. S12 we show the effect of only including the inner four rings (1 • < r < 4 • ) in our analysis for the NFW profile (left) and Romulus profile (right). The results are consistent with those found in Fig. 3, as would be expected considering Fig. S9, though the best-fit cross-sections are slightly larger than in the analyses that include all annuli.
In Fig. S13, on the other hand, we remove the inner four rings and only keep the outer five rings, such that 4 • < r < 10 • . The results are consistent with those in Fig. 3 and the inner Galaxy results in Fig. S12, with the best-fit cross-sections being slightly smaller.
In Fig. S14 we increase the plane mask to only include |b| ≥ 2 • . As the result of the increased plane mask the innermost ring is completely masked, and so we only include annuli 2 through 9. This analysis is meant to address possible systematic effects from mismodeling near the Galactic plane, though the results are consistent with those of the fiducial analysis.

B. Energy range variations
We now revert to our fiducial analysis ROI but we change the number of energy bins included in the spectral fits. In particular, we consider adding in one additional (lower) energy bin, so that photons with E 8 GeV are included, and then removing 3 energy bins so that only photons with E 20 GeV are incorporated. Then, we remove highenergy photons so that the energy constraint is 10 GeV < E < 100 GeV. The results of these analysis variations are illustrated in Fig. S15. Including lower-energy photons produces consistent results to our fiducial analysis but slightly increases the discovery TS, while removing the first few energy bins (E > 20 GeV) lowers the TS but still produces a consistent best-fit and upper limit. Removing the upper energy bins (E < 100 GeV) also leads to consistent results, as illustrated in the bottom panel.
As the lower energy cut-off is increased, we rapidly lose sensitivity to a putative higgsino signal. As an illustration, in Fig. S16 we show the results of analyses where we restrict E > 40 GeV and E > 50 GeV. Interestingly, there is still a ∼1σ preference for a higgsino, though the uncertainties are significantly increased compared to in our fiducial analyses. . As in Fig. 3 and Fig. S12 but with the Galactic plane masked at |b| ≥ 2 • . The innermost ring is not included since it is fully masked.

C. Removing the PS and isotropic spectral templates
In our fiducial analysis we include Galactic emission, PS, and isotropic spectral templates. However, as shown in Fig. S1, the PS and isotropic templates are subdominant compared to the Galactic emission in our energy range and ROI. As a cross-check of our analysis, we consider the effects of removing the PS and isotropic templates completely from our analysis. The results of this analysis are shown in Fig. S17. Comparing these results to Fig. 3, the differences are minor and within statistical uncertainties.

D. Data-driven background model
In our fiducial analysis we use the p8r3 Galactic emission model in addition to the PS and isotropic models to describe the null hypothesis. In this section we instead adopt a data-driven approach and determine the background model from the ROI described by 10 • < r < 15 • , with PSs masked (using our fiducial mask) and the Galactic plane masked such that |b| ≥ 1 • . In the background ROI we determine the flux spectrum dF/dE numerically, in units of cts/cm 2 /s/sr/GeV, and we assume the background emission within our signal ROI takes on the same spectral shape. However, we still assign a nuisance parameters that re-scale the overall background flux dF/dE in each annulus independently. Crucially, we do not make use of the p8r3 Galactic emission model in this analysis, and we make no assumptions about the spatial dependence of the background model. We also do not include separate PS and isotropic spectral emission templates, since these should be captured by the data-driven template.
Given that PSs and Galactic emission are more important at low energies, we restrict this analysis to E > 20 GeV. We also restrict E < 200 GeV because of the large statistical uncertainties in determining the background-region model at high energies. The best-fit data-driven model is shown relative to the signal region data (summed over all 9 annuli) in Fig. S18. In that figure we also compare the spectrum to our best-fit fiducial model (labeled p8r3) which is also determined in our fiducial energy range. Note that while in this figure we show the summed spectrum over the full ROI the analysis is performed through a joint likelihood over the independent annuli, as in our fiducial analysis.
We do not account for statistical uncertainties in our determination of the background-region spectrum in our analysis, since the background region is larger than the signal region. The statistical uncertainties on the backgroundregion spectrum are illustrated in Fig. S18 relative to the summed signal region statistical uncertainties. As a further cross-check, we consider only including the inner four annuli (1 • < r < 5 • ), such that the signal region is much smaller than the background region and such that the two regions are separated by 5 • . Including all annuli out to 10 • , the best-fit cross-section for the NFW (Romulus) profile at the thermal mass is 10 ± 9 × 10 −26 cm 3 /s (2.2 ± 1.5 × 10 −26 cm 3 /s), while restricting to the inner four annuli the best fit changes to 18±12×10 −26 cm 3 /s (3.0±1.8×10 −26 cm 3 /s). Thus, shrinking the signal region for the data-driven analysis leads to consistent results as in the full signal-region analysis. In Fig. S19 we show the results of the analysis searching for the higgsino signal (including all signal-region annuli) with the data-driven background model. The results are consistent with those in Fig. 3, though the search is somewhat less constraining because of the reduced energy range.

E. Including all data quartiles
In this section we increase the data volume relative to our fiducial analysis. In the main Letter we use the top 3 quartiles of events as ranked by PSF. Here, we include all 4 quartiles of events. Naively including more events should increase the sensitivity to a higgsino signal (with e.g. detection significances increasing, on average, by around 15% in the event of a signal), however this is less straightforward in our analysis framework since including the last quartile of events increases the size of the PSF mask.
In Fig. S20 we show the results of the higgsino search using the full data set, as presented in Fig. 3 for our fiducial analysis. The results including the last quartile of data are consistent with those in our fiducial analysis.  Figure S18. The best-fit spectra for our fiducial analysis and the data-driven background model analysis, which is described in Sec. III D. Note that we restrict the energy range to 20 GeV < E < 200 GeV for the data-driven background model analysis. The spectra are shown summed over all 9 annuli, over the full ROI, though the analysis uses the joint likelihood constructed by analyzing the annuli independently. We indicate statistical uncertainties in our determination of the background model, which are subdominant compared to the statistical uncertainties in our summed signal region by a factor ∼0.8 across the energy range.  Fig. 3 and Fig. S12 but with a data-driven background model that consists of the spectrum from 15 • > r > 10 • with |b| < 1 • masked in addition to PSs masked. We restrict the energy range to 20 GeV < E < 200 GeV for this analysis.

F. Including an inverse Compton template
In our fiducial analysis we use the p8r3 Galactic emission model, which is constructed from multiple components. Broadly speaking, p8r3 includes gas-correlation emission templates in addition to IC templates that correlate with the radiation field and the cosmic ray distribution. Since these two sources of emission have relatively uncorrelated systematics, it makes sense to try to assign them individual nuisance parameters in the analysis. We implement this by using the p8r3 emission template in addition to an extra IC spectral template, which is allowed to have an unconstrained positive or negative normalization. We use the IC template from Model F of [60], reprocessed for our data set. The IC template is illustrated in Fig. S1. The results of including this template are shown in Fig. S21. The best-fit cross-section is consistent with our fiducial result, though slightly larger.  Fig. 3 and Fig. S12 but including an additional IC template that may take on a positive or negative normalization.

G. Spatial template fit analysis
We now present the results of the spatial template fit analysis, which is described in Sec. I and which makes use of the likelihood in (S3). As a reminder, we include three spatial templates in each energy bin over our ROI, with independent nuisance parameters in each energy bin that rescale the overall template normalizations. The templates are Galactic emission (p8r3), PS emission from the 4FGL that extends beyond our PS mask, and DM annihilation as given by the J-factor map. Our analysis ROI is a slight variation of that found by joining all of our annuli, as described below. We consider all energy bins above 10 GeV and below the DM mass, as in our fiducial spectral analysis.
The spatial template analysis is more susceptible to mismodeling since we now have to describe the data in each spatial pixel in addition to each energy bin. The Galactic emission templates are notoriously poor descriptions of gamma-ray emission in the inner Galaxy (see, e.g., [61]). Indeed, we find (by visual inspection) that the analysis in our fiducial ROI produces a poor fit to the data, with large residuals and regions of over-subtraction near the Galactic plane in particular. We partially mitigate the mismodeling by increasing the size of our Galactic plane mask, increasing it to |b| ≥ 1.5 • and to |b| ≥ 2 • . Of course, in increasing the plane mask we also mask more region near the GC, which decreases our sensitivity to a putative DM signal. In Fig. S22 we show the best-fit null-hypothesis models (left panels) and smoothed residuals (right panels) from analyses with |b| ≥ 1.5 • (top panels) and |b| ≥ 2 • (bottom panels). While we perform the analyses in each energy bin independently, for these figures we sum the results over all energy bins (taking m χ = 1.1 TeV so that we only include photons with E 1.1 TeV).
In the top right panel of Fig. S22 large regions of over-subtraction are visible by eye, particularly in the northern hemisphere near the Galactic plane. The over-subtraction is better though still present in the |b| ≥ 2 • analysis.
The results of the spatial template analysis (to be compared to Fig. 3) are illustrated in Fig. S23. Interestingly, the results of this analysis are broadly consistent with those found in the fiducial spectral analysis, with comparable upper limits and consistent (positive) best-fit values, though the detection significance in slightly reduced in the spatial analyses. The results of the |b| ≥ 1.5 • and |b| ≥ 2 • analyses are consistent, though as expected the |b| ≥ 1.5 • analysis is slightly more sensitive. Note that we only present results using the NFW DM profile since the FIRE-2 profiles are not at high enough spatial resolution to be used in the template analysis without further smoothing.   Figure S23. As in Fig. 3 and Fig. S12 but for the spatial template analyses using the NFW DM profile with |b| ≥ 1.5 • and |b| ≥ 2 • , as indicated.