Faint Dark Matter Annihilation Signals and the Milky Way's Supermassive Black Hole

A wide range of mechanisms predict present-day s-wave dark matter (DM) annihilation cross-sections that are orders of magnitude below current experimental sensitivity. We explore the capability of DM density spikes around the Milky Way's supermassive black hole to probe such faint signals of DM annihilations, considering a range of possible spike and halo distributions. As an exemplar of a theory with a suppressed s-wave annihilation cross-section, we consider a hidden sector axion portal model of DM. In this model, the leading contribution to the annihilation cross-section in the early universe is p-wave, while s-wave annihilations occur at higher order in the coupling constant. We provide a unified treatment of DM freezeout in this model including both s- and p-wave annihilations and analytically determine the photon spectrum for the dominant DM annihilation process in the universe today. We find that Fermi and H.E.S.S. observations of the Galactic Center offer excellent sensitivity to this model over a wide range of parameter space, with prospects depending sensitively on the properties of the DM spike as well as the central halo.


I. INTRODUCTION
The indirect detection of dark matter through its annihilation products in cosmic rays is a cornerstone of the experimental search for dark matter (DM). Indirect detection is an increasingly potent probe of annihilating DM, with observations of (e.g.) both the cosmic microwave background and dwarf galaxies now sensitive to DM with annihilation rates at or below the standard thermal target σ thermal = 3 × 10 −26 cm 3 /s in an expanding range of masses and annihilation channels [1][2][3][4].
Identifying the annihilation products of TeV-scale DM with standard thermal-scale cross-sections remains a steep observational challenge, however, as the flux of cosmic rays from DM annihilation in galaxy halos falls off with increasing DM mass as m −2 . Moreover, many dark matter models predict present-day DM annihilation cross-sections substantially below the thermal target. There are many mechanisms that predict a suppressed present-day DM annihilation cross-section. For instance, the DM may simply arise from a dark sector that is very cold compared to the Standard Model (SM) [5,6]; it may have been diluted by entropy production post-freezeout [7][8][9]; its annihilation cross-section may be strongly velocity-dependent, because of kinematics [10][11][12][13] or symmetries [14][15][16][17][18] (or both); or its annihilation cross-section in the early universe may involve additional species that are later depleted [10,[19][20][21][22]. In several of these models, e.g. [5,6,17], terrestrial signals are typically significantly reduced compared to expectations from thermal WIMPs, making even a suppressed indirect detection signal an irreplaceable discovery handle and a powerful window onto the physics of DM.
Here we estimate the sensitivity to faint s-wave DM annihilation cross-sections that can potentially be offered * tchiang5@illinois.edu † slshapir@illinois.edu ‡ sheltonj@illinois.edu by DM density spikes around the supermassive black hole (SMBH) at the center of the Milky Way. Black holes focus DM within their gravitational zone of influence into a steep, localized over-density known as a spike [23]. This enhancement of the DM density can potentially magnify DM annihilation rates by many orders of magnitude in the immediate vicinity of the black hole, leading to a bright, point-like source of cosmic rays. While this enormous magnification of DM annihilation signals can offer a uniquely powerful window onto models of DM with suppressed annihilation cross-sections and thus suppressed signatures in DM haloes [17,[24][25][26][27], the details of the predicted spike distribution depend sensitively on as-yetunknown properties of the host DM halo, the central stellar cusp within that halo, and the formation history of the black hole [23,[28][29][30][31][32][33]. In this work we consider a broad range of possible DM distributions in the Galactic Center, with the aim of understanding what parameter ranges offer interesting sensitivity to sub-thermal annihilation cross-sections.
As a representative DM model with challenging annihilation cross-sections, we consider the hidden sector axion portal (HSAP) model developed in [17,27]. In this model, fermionic DM χ annihilates to pseudo-scalars a, which subsequently decay to SM electroweak gauge bosons. This model features an interesting interplay of two annihilation channels: the process χχ → aa proceeds in the p-wave, while the reaction χχ → aaa contributes in the s-wave when it is kinematically available, but is higher order in the coupling constant. For low DM masses, the p-wave process dominates DM annihilation in the early universe, yielding a subdominant s-wave annihilation cross-section orders of magnitude below the standard thermal target. For high DM masses, the s-wave annihilation process can dominate during thermal freezeout, resulting in thermal-scale cross-sections but signals that are observationally challenging thanks to the large DM mass. In all cases the s-wave process χχ → aaa dominates the present-day DM annihilation rate, including within DM spikes.
We begin in Section II with an overview of the HSAP model, including a novel treatment of freezeout incorporating both s-and p-wave contributions to the annihilation cross-section, and a calculation of the photon spectrum from the resulting DM annihilations. In Sec. III we detail our model of DM density spikes around the Milky Way's SMBH. We compare predicted gamma-ray fluxes from DM annihilation within SMBH-induced density spikes to observations from Fermi and H.E.S.S. in Sec. IV and discuss the resulting prospects for sensitivity, and in Sec. V we conclude.

II. THE PARTICLE MODEL
In this section we discuss a hidden sector axion portal model of dark matter, as introduced in [17]. In this model, DM is a Majorana fermion χ which annihilates to pseudoscalars a, which subsequently decay to the SM via axion-like couplings to SM gauge bosons. This model is CP -conserving, ensuring that the leading annihilation process χχ → aa is p-wave. In contrast to previous works [17,27], we will focus on the regime where the higherorder but s-wave annihilation process χχ → aaa is kinematically available, and explore the consequences both for thermal freezeout and for potential BH spike signals in our Galaxy today.

A. Annihilations and Relic Abundance
The hidden sector axion portal (HSAP) model is described by the Lagrangian This Lagrangian has three free parameters: the masses m χ and m a , and the Yukawa coupling constant y, which can be determined in terms of m χ and m a using the requirement that thermal freezeout of DM annihilations yields the observed DM relic abundance. Additionally, the mediator a is coupled to the SM via dimension-five axion-like interactions with SM gauge fields, which enables it to decay promptly on astrophysical scales, as we discuss further below.

Annihilation cross-sections
The leading 2 → 2 annihilation process occurs in the p-wave, with thermally averaged cross-section given in the non-relativistic limit by [17] σv p = 1 x where η ≡ m a /m χ and x ≡ m χ /T.
The velocity dispersion is related to the temperature as v 2 = 6/x.
There is one Feynman diagram topology for the process χχ → aaa, shown in Fig. 1, which yields six distinct diagrams thanks to the permutation of final state momenta among the identical particles. The calculation of the matrix element simplifies significantly in the nonrelativistic limit, where all of the angular dependence drops out. The resulting expression for the spin-averaged matrix element is where the variables w i parametrize the distribution of energy among the three final state particles in the CM frame, and (7) In the non-relativistic limit, the cross-section for χχ → aaa is then where the upper and lower limits of integration for w 2 are given by and for w 1 the upper and lower limits are The velocity-dependent O(v 2 ) term is negligible both at freezeout (where it is higher order in y 2 /4π compared to the contribution from χχ → aa) and in the Milky Way today. We thus retain only the piece of Eq. 8 that is constant as v 2 → 0, defining a contribution to the s-wave annihilation σv s .

Thermal freezeout and relic abundance
To determine the Yukawa coupling constant y, we include the leading contributions to both s-and p-wave annihilation processes and solve the Boltzmann equation governing the DM relic abundance. In the nonrelativistic limit, this Boltzmann equation can be written as where λ = 0.264(g * S /g 1/2 * )m Pl m χ and the equilibrium yield is Y eq = 0.29(g * S ) −1 x 3/2 e −x [34]. Here m Pl = G −1/2 is the Planck mass, and g * and g * S are the number of effective relativistic degrees of freedom contributing to the energy and entropy densities, respectively. The present-day DM relic abundance is taken to be Ω DM h 2 = 0.112 [1]. Figure 2 shows the resulting contours of y 2 /4π in the m χ -η parameter space. The model becomes nonperturbative for y 2 /4π 1, which restricts m χ 14 TeV. The interplay between s-and p-wave contributions during freezeout becomes especially important for heavier m χ and smaller η, where the χχ → aaa annihilation is less suppressed compared to the χχ → aa process. Figure 3 compares the s-wave (solid) and p-wave (dashed) annihilation cross-sections for four representative values of the mass ratio, η = 0.05 (blue), 0.3 (yellow), 0.6 (green), and 0.66 (red). In this figure we evaluate the p-wave annihilation cross-section at freezeout, x f , as determined through the sudden freezeout condition Here the Hubble rate is given by H(x f ) = The resulting values of x f range from 22 < x f < 32 for 5 GeV < m χ < 14 TeV and 0 < η < 2/3. Since x f increases with increasing m χ , at large m χ the s-wave cross-section becomes increasingly important at x f compared to the velocity-suppressed pwave contribution, and can dominate freezeout for sufficiently heavy DM and sufficiently small η. The size of the s-wave cross-section depends sensitively on the mass ratio η, and is strongly suppressed as η approaches the kinematic limit of 2/3.
Since the typical DM velocity dispersion in the Galaxy today is v 2 = O(10 −6 ), whereas v 2 = O(10 −1 ) at x f , the p-wave cross-section today is suppressed by five to six orders of magnitude relative to its value at freezeout. Figure 3 thus demonstrates that DM annihilations in the galaxy today are dominated by the s-wave χχ → aaa process. We observe two distinct regimes, depending on whether the p-or s-wave process dominates at x f . When the p-wave process dominates, at small DM mass and large η, the s-wave annihilation cross-section is suppressed by orders of magnitude compared to the typical thermal target (∼ 1 pb). On the other hand, when the s-wave cross-section dominates, at large DM mass and small η, it is comparable to the thermal target.

Pseudoscalar decays into SM final states
The pseudoscalar a can decay to the SM through dimension-five axion-like couplings to SM gauge bosons. For simplicity, we will consider here the case when the pseudoscalar couples at leading order only to the hypercharge field strength B µν via the interaction This choice is also a relatively optimistic scenario for discoverability, as it yields an energetic spectrum with a distinctive spectral feature. After electroweak symmetry breaking, the interaction of Eq. 13 mediates the decays a → γγ and, if kinematically possible, a → Zγ and a → ZZ. For a given value of m a , the partial widths into γγ, Zγ and ZZ final states following from the interaction of Eq. 13 are given by where θ W is the Weinberg angle and C ∼ m a /Λ 2 is a common constant of proportionality. We require C to be small enough that annihilations of DM directly to SM gauge bosons through an intermediate a are negligible in comparison to the secluded annihilation processes of Fig. 1, but otherwise our results do not depend on C.
We require only the relative branching fractions of the pseudoscalar into the various decay channels that follow from Eq. 14.

B. Photon Spectra
In this subsection we determine the photon spectra dN/dE γ resulting from DM annihilations χχ → aaa.

Photon spectra from a → γγ decays
We begin with the photon spectrum dN/dE γ from an s-wave annihilation process χχ → aaa → 6γ, which can be obtained analytically in the non-relativistic limit. Consider an individual a particle with energy E a in the As η increases toward the kinematic limit, the restricted phase space forces the spectrum to become increasingly narrow and peaked.
Galactic frame decaying into two photons. The maximum and minimum energies the photons can have are As the decays of a are isotropic in its rest frame, the energy distribution of daughter photons is uniform between the kinematic boundaries. Defining the dimensionless variables and w ≡ E a /m χ , the kinematic endpoints can be written For a given w and η, the probability P (u; w, η)du of finding a photon within the energy interval du is thus as the probability is uniform and unit-normalized over the kinematically allowed interval. The probability of obtaining an a particle with (relative) energy w from the process χχ → aaa is given by where the cross-section σ for χχ → aaa is given in Eq. 8 and A is a normalization factor. The probability density of finding any specific combination of u, w is then As this expression indicates, a daughter photon with energy u may have come from a parent a with a range of possible energies w. To obtain the probability of observing a photon with energy u, we integrate P (u, w; η) over the range of w consistent with the value of u. The maximum possible value of w is, from Eq. 10, w max = 1 − 3 4 η 2 , independent of u. For a fixed value of u, w min can be determined from Eq. 17, which gives The desired photon spectrum is therefore The factor of six appears here since there are six final state photons in the annihilation χχ → aaa → 6γ. Accordingly, we evaluate the normalization factor A by requiring In Fig. 4, the normalized photon spectra, weighted by u 2 , are plotted for four different values of η. For small η, the energy distribution is broad and peaks at relatively high energies compared to the spectra for larger values of η. As η approaches the kinematic limit of 2/3, the changes in the spectrum shape become increasingly rapid as the available phase space shrinks.

Photon spectra from a → Zγ and a → ZZ decays
Once m a > m Z (2m Z ), the decay channel a → Zγ (a → ZZ) opens up. These decay channels result in a continuum of numerous but lower-energy photons from hadron decays, as well as final state radiation off of charged leptons. To obtain photon spectra dN/dE for DM annihilation channels with any number of finalstate Z bosons, we proceed numerically. A Monte Carlo sampling method was employed to compute the fourmomenta of a particles produced in χχ → aaa annihilations following the non-relativistic distribution of Eq. 7, as well as the momenta of their daughter photons and Z bosons. The photon spectrum resulting from a Z boson in its rest frame was computed using Pythia 8 [35]. The resulting photons were then boosted to the Galactic rest frame. Figure 5 shows the normalized photon spectra resulting from a decays to γγ (yellow), ZZ (blue), and Zγ (green) for an example parameter point with m χ = 1500 GeV and η = 0.3. The γγ decay process produces comparatively more high-energy photons with a narrower distribution. The spectrum from ZZ decay is considerably broader with a smooth low-energy tail. The Zγ decay spectrum inherits features from both γγ and ZZ spectra, albeit with a slightly smaller maximum photon energy than the γγ spectrum due to the non-zero m Z .

III. THE ASTROPHYSICAL MODEL
Following [17,27,36], we consider a simple parametric model describing possible DM density spikes at the Galactic Center. We take the DM halo to be described by a generalized Navarro-Frenk-White (NFW) model, which in the central regions of the Galaxy, i.e., well within the scale radius, can be described by a power law ρ(r) = ρ(r 0 )(r/r 0 ) γc . We anchor this power law to the Solar System, where we take the local density of DM to be ρ(r ) = 0.3 GeV/cm 3 [37], at a distance r = 8.46 kpc from the Galactic Center [38]. DM-only simulations typically yield values of the cusp exponent γ c within the range 0.9 γ c 1.2 [39,40], while the adiabatic contraction of the central halo following baryonic collapse can lead to larger values of γ c [41][42][43]. We will treat γ c as a free parameter.
The Milky Way's SMBH has mass M = 4 × 10 6 M [44,45]. Its gravitational zone of influence extends out to the radius r h ≡ GM/v 2 0 , where the gravitational potential energy due to the BH is equal to the typical kinetic energy of a DM particle in the halo. Here v 0 is the velocity dispersion in the inner halo. We adopt as our fiducial dispersion v 0 = 105 ± 20 km s −1 [46]. Numerical studies indicate that the spike itself begins growing somewhat within r h , at r b = 0.2r h [30,47]. The spike is also well-described as a power law, ρ sp (r) = ρ(r b )(r/r b ) γsp , although the spike index γ sp depends sensitively on the formation history of the SMBH and the properties of its environment. Spikes that form around a BH that grows adiabatically at the center of a cuspy DM halo are very steep, γ sp = (9 − 2γ c )/(4 − γ c ) [23]; on the other hand, if the BH is not at the dynamical center of its halo, then it produces a very shallow spike, γ sp = 1/2 [28,29]. The dynamical heating of DM from gravitational interactions with a dense and cuspy stellar distribution results in a spike solution with limiting index γ sp = 1.5, attained when the system has fully equilibrated [30][31][32]. Non-equilibrated spikes with intermediate power laws are possible if the DM at the Galactic Center is still in the process of equilibrating. 1 Here we treat γ sp as a free parameter and vary it between 1.5 and its adiabatic value.
Once the DM density in the spike reaches the value ρ ann = m χ /( σv τ ), where τ ≈ 10 10 yr is the lifetime of the spike, DM annihilations become rapid enough to deplete the spike. We define the radius at which this occurs as r in = r b · (ρ b /ρ ann ) 1/γsp . For r < r in , the spike follows a very mild power law, ρ(r) ∝ r −1/2 in the case of s-wave annihilations [50,51]. Finally the inner boundary of the spike is obtained at r < 4GM [52]. Altogether, then, we take for our model where ρ sp (r) = ρ b (r b /r) γsp , ρ in (r) = ρ ann (r in /r) 0.5 , ρ b = ρ(r ) · (r /r b ) γc , and we we have defined a spike profile that smoothly interpolates between the inner spike with index 1/2 and the outer spike with index γ sp . With the parameter values adopted here, r b ≈ 0.3 pc, while for typical DM parameters in the HSAP model the radius r in ∼ few × 10 −5 pc. The inner spike structure, which dominates the emission, then subtends ∼ 4 milliarcseconds on the sky, which is several orders of magnitude below the typical resolution of the Fermi telescope. We thus treat the spike signal as a point source.
Given a photon spectrum dN/dE, the differential photon flux observed on Earth is then [53] where the line-of-sight (l.o.s.) distance at an angle θ away from the GC is given by We conservatively include the emission from the central 0.5 • of the halo, and indicate where this halo emission, Φ halo , exceeds the contribution from the spike, Φ sp . Our 1 An alternate parameterization of non-equilibriated spikes undergoing a baryonic heating process, developed in [48], gives indistinguishable results for the spike signal [49]. main interest is in the regime where the spike dominates, Φ sp Φ halo . This choice allows us to succinctly account for both the finite angular resolution of gamma-ray telescopes, notably Fermi [54], and demarcate the region where we might expect challenges from resolving a point source on top of a bright and non-uniform background.
The spike component is treated as a point source, as most emission comes from the innermost region within radius r in in this model. The spike contribution is then given by

IV. SENSITIVITY TO FAINT s-WAVE ANNIHILATION SIGNALS
Both Fermi and H.E.S.S. have observed bright point sources in the Galactic Center that the respective collaborations have associated with Sgr A*. The Fermi point source 3FGL J1745.6-2859c is observed in gamma rays with 100 MeV < E γ < 100 GeV [55], while the H.E.S.S. source HESS J1745-290 is observed in gamma rays with 180 GeV < E γ < 79 TeV [56]. We show the observed spectra of these point sources in the solid histograms of Figure 6. The abrupt decrease in binned flux magnitude starting at m χ = 180 GeV reflects the different bin sizes in the two experiments. The energy gap between Fermi and H.E.S.S. datasets at 100 < m χ < 180 GeV is clearly visible. For comparison, the dashed histograms show four example predictions for the primary photon spectra arising from HSAP DM annihilations within DM density spikes. We show predictions for four different values of DM mass at fixed η = 0.1, γ c = 1.2, and γ sp = 1.8.
To set constraints on HSAP DM annihilations within BH spikes, we use the simple criterion that the spike flux does not exceed the observed flux from either 3FGL J1745.6-2859c or HESS J1745-290 at more than 95% CL in any bin. Figure 7 shows the values of the s-wave annihilation cross-section excluded by this procedure for six fixed values of m χ in two different astrophysical scenarios. In this figure we take Br(a → γγ) = 1 throughout, so that the photon spectrum is given by Eq. 23. Visible kinks in the lines reflect when the peak of the photon spectrum moves from one bin to another as η changes. For m χ = 5, 11, and 50 GeV, the constraints come from Fermi, while for m χ = 750, 1500, and 6500 GeV, the constraints come from H.E.S.S. The improvement of the limits as η → 2/3 reflects the increasing sharpness of the peak in the photon spectrum. From Fig. 7, we see that adiabatic spikes are sensitive to cross-sections some four orders of magnitude smaller than the standard thermal target, while for more moderate values of γ sp , the resulting sensitivity even in relatively cuspy haloes is less dramatic but still interesting. Constraints on the crosssection for higher mass DM are generally weaker; however, in the HSAP model studied in this paper, this effect is substantially offset by the larger predicted s-wave cross-sections at high mass. Figure 8 shows excluded regions of m χ -η parameter space for fixed choices of spike and halo indices. The top panel fixes γ c = 1.1 together with γ sp = 1.9, 2.0 (no parameter space is excluded for γ sp ≤ 1.8). The bottom panel has γ c = 1.2 with γ sp = 1.7, 1.8, 1.9, 2.0. For each exclusion boundary, a concave curve segment corresponds to the exclusion given by one energy bin from either the Fermi or H.E.S.S. point source. As expected from Fig. 3, we see strong variation of the sensitivity with the mass ratio η at small m χ .
In Fig. 8 we continue to take for simplicity Br(a → γγ) = 1, as this case can be handled analytically. We indicate with the shaded light (dark) tan area where m a > m Z (2m Z ), and this assumption breaks down. In general, once Z bosons are kinematically accessible, the details of the DM annihilation spectrum will depend on the modeling choice of how the pseudoscalar a couples to SM electroweak bosons. However it is straightforward to translate the results of Fig. 8 to the non-zero branching fractions into Zγ and ZZ predicted by Eq. 14 in most of parameter space, as we now discuss.
Exclusions in the HSAP model are dominated by the peak of the photon spectrum in the γγ and Zγ channels; the continuum photons coming from Z bosons are not important for determining the sensitivity. This can be seen from Fig. 9, where we compare limits set using only γγ decays to those determined from both γγ and Zγ decays. This figure shows the maximum s-wave annihilation cross-section allowed by observations of the GC, relative to the prediction of the HSAP model, as a function of DM mass, and demonstrates how the limit changes as the branching fractions of a are varied. Green, yellow, and blue curves show exclusions assuming branching ra- In the bottom panel, the blue curves use the branching ratios of Eq. 14, while the green curves, marked with asterisks, artificially assume Br(a → γγ) = 1.
i.e., by using the γγ limit but multiplying by the factor Br(γγ) + 1 2 Br(Zγ) < 1. This is demonstrated in Fig. 9 by the excellent agreement of the blue dashed and orange solid curves at large m χ . Figure 9 also shows the expected linear scaling of the high-m χ excluded cross section with Br(γγ). By contrast, in the gap between Fermi and H.E.S.S.' observations, the limit is dominated by the high-energy shoulder of the photon spectrum, where photons from Zγ decays make a negligible contribution.
As Fig. 8 already demonstrates, observational sensitivity to the HSAP model is strongly dependent on γ sp and especially γ c . Any astrophysical parameter combination with γ c ≥ 1.2 and γ sp ≥ 2.0 is almost entirely ruled out. Conversely, very little parameter space is excluded for γ c ≤ 1.1 and γ sp ≤ 1.8. We now turn to examining the sensitivity as a function of γ c and γ sp for fixed particle parameters in Fig. 10. Here we show results for fixed m χ = 11 GeV (top) and 1500 GeV (bottom), and three fixed choices of mass ratio η. The gray shaded regions indicate where the contribution from the central 0.5 • of the halo exceeds the contribution from the spike. The exclusions we find are thus indeed driven by the spike for the bulk of parameter space, except at large values of γ c where a more careful combined study of the halo plus potential spike point source would be warranted. Canonical adiabatic spikes are almost entirely disfavored, except at very small values of γ c .
In general, points with η = 0.1 see more stringent limits than those with η = 0.3 and 0.6. This occurs for two main reasons. First, in the HSAP model, the magnitude of the s-wave cross-section depends sensitively on η (see Fig. 3), especially for the relatively low value m χ = 11 GeV in the top panel. Second, once it is kinematically possible to produce Z bosons, the flux in the photon peak is reduced. The branching fraction γγ+ 1 2 Zγ is typically larger for smaller values of η. This second effect is relevant in the bottom panel of Fig. 10, where the pseudoscalar is heavy enough to decay into Zγ (η = 0.1, 0.3) and ZZ (η = 0.6) final states. In this panel we show two sets of limits. The green curves, marked with asterisks, show limits obtained by artificially setting Br(γγ) = 1, while the blue curves show limits using the branching fractions of Eq. 14. In particular, we set {Br(a → γγ), Br(a → Zγ)} ≈ {0.65, 0.35} for η = 0.1, {Br(a → γγ), Br(a → Zγ), Br(a → ZZ)} ≈ {0.32, 0.61, 0.07} for η = 0.3, and {Br(a → γγ), Br(a → Zγ), Br(a → ZZ)} ≈ {0.3, 0.62, 0.08} for η = 0.6. Single asterisks indicate that a → Zγ is kinematically allowed, while double asterisks indicate both a → Zγ and a → ZZ decay channels are possible.
The exclusions for η = 0.1 and 0.3 in the bottom panel of Fig. 10 are set by the same bin, while for η = 0.6 the peak of the photon spectrum has migrated far enough to higher energies that it is instead the neighboring bin that determines the exclusion. This results in the different shapes of the excluded regions obtained for η = 0.6 versus η = 0.1 and 0.3 that can be most easily seen in the green curves. When we incorporate the above branching ratios into Zγ and ZZ, the flux in the peak of the photon spectrum is reduced and the limits weaken accordingly, but the shapes of the curves remain the same.

V. SUMMARY AND CONCLUSIONS
We have performed a detailed study of DM annihilations in a hidden sector axion portal (HSAP) model where fermionic DM χ annihilates to pseudoscalars a. We consider here the regime where both the p-wave χχ → aa and s-wave χχ → aaa annihilation channels are kinematically available. We find that the relic abundance is determined dominantly by p-wave annihilations at low DM masses, while s-wave contributions can dominate at higher DM masses. In all cases the s-wave crosssection dominates the annihilation in the Milky Way today, and can be suppressed by as much as 10 −4 relative to the standard expectation for a thermal WIMP. We take the pseudoscalar mediator a to decay to SM states via an axion-like coupling to SM hypercharge gauge bosons. When m a < m Z , the pseudoscalar decays exclusively via a → γγ. We derive the photon spectrum resulting from χχ → aaa → 6γ analytically in the limit of nonrelativistic DM annihilations. This spectrum has a characteristic high-energy spectral feature that could help to facilitate discovery and identification of a DM signal in the busy environment of the Galactic Center. When m a is large enough to allow for Zγ and ZZ decays, we determine the spectra for those decay channels numerically.
DM density spikes around the Milky Way's supermassive black hole offer a speculative but irreplaceable potential discovery handle on this model. We estimated current sensitivity to HSAP DM annihilation within such spikes by comparing the predicted point-like gamma-ray spike signals to Fermi and H.E.S.S. observations of Sgr A*. We find interesting sensitivity to the HSAP model across a wide range of both particle and astrophysical parameters. Notably, canonical adiabatic spikes can probe cross-sections some four orders of magnitude below the standard thermal target, and are accordingly almost entirely ruled out in our model space. However substantial portions of parameter space remains open for shallower spikes, particularly in less cuspy haloes.
The HSAP model compactly provides examples of both parametrically suppressed annihilation cross-sections at low DM mass and mass-suppressed annihilation rates at high DM mass. As is common for models of secluded DM [57], the HSAP model also has severely suppressed signals at both direct detection and collider experiments, making indirect detection a vital experimental probe of this model. This HSAP model is only one example of a rich variety of particle models that predict DM annihilation signals substantially below current observational sensitivity; the enhanced sensitivity to such faint annihilation signals in DM density spikes thus remains important to consider even as thermal targets are surpassed.