Prospects for indirect dark matter searches with MeV photons

TL;DR

This paper evaluates the potential for future MeV gamma-ray missions to detect dark matter annihilation signals in the 1 MeV to 3 GeV mass range, emphasizing spectral features like bremsstrahlung and in-flight annihilation.

Contribution

It provides a comprehensive analysis of detection prospects for sub-GeV dark matter in the MeV range, focusing on spectral signatures often neglected in previous studies.

Findings

Bremsstrahlung and in-flight annihilation are key spectral features for detection.

Future MeV missions could significantly improve sensitivity to sub-GeV dark matter.

Leptonic channels offer promising detection signatures.

Abstract

Over the past decade, extensive studies have been undertaken to search for photon signals from dark matter annihilation or decay for dark matter particle masses above $\sim1$ GeV. However, due to the lacking sensitivity of current experiments at MeV-GeV energies, sometimes dubbed the 'MeV gap', dark matter models with MeV to sub-GeV particle masses have received little attention so far. Various proposed MeV missions (like, e.g., e-ASTROGAM or AMEGO) are aimed at closing this gap in the mid- or long-term future. This, and the absence of clear dark matter signals in the GeV-TeV range, makes it relevant to carefully reconsider the expected experimental instrumental sensitivities in this mass range. The most common two-body annihilation channels for sub-GeV dark matter are to neutrinos, electrons, pions or directly to photons. Among these, only the electron channel has been extensively studied, and almost exclusively in the context of the 511 keV line. In this work, we study the prospects for detecting MeV dark matter annihilation in general in future MeV missions, using e-ASTROGAM as reference, and focusing on dark matter masses in the range 1 MeV-3 GeV. In the case of leptonic annihilation, we emphasise the importance of the often overlooked bremsstrahlung and in-flight annihilation spectral features, which in many cases provide the dominant gamma-ray signal in this regime.

Figures

Tables

Table 1. Table 1 : Instrumental details for three proposed γ 𝛾 \gamma –ray telescopes, e-ASTROGAM [ 10 ] , COMPAIR [ 11 ] and ADEPT [ 12 ] . In this work we assume an on-target observation time of T obs = 1 ​ yr subscript 𝑇 obs 1 yr T_{\mathrm{obs}}=1\mathrm{\,yr} . Our reference experimental scenario is e-ASTROGAM performing a full-sky survey. All details are taken from Ref.  [ 10 ] . The energy resolution in the pair conversion regime is set to Δ ​ E / E = 30 % Δ 𝐸 𝐸 percent 30 \Delta E/E=30\% .

Compton domain
	e-ASTROGAM	COMPAIR	ADEPT
Energy range $\left[\mathrm{MeV}\right]$	0.3–10	0.2 – 10	–
$\Delta E/E$	1.3%	2%–5%	–
$A_{\mathrm{eff}}\left[{\mathrm{cm}}^2\right]$	50–560	50 – 250	–
FoV $\left[\mathrm{sr}\right]$	2.9	3	–
Pair-conversion domain
	e-ASTROGAM	COMPAIR	ADEPT
Energy range $\left[\mathrm{MeV}\right]$	10–3000	10–500	5–200
$\Delta E/E$	20–30%	12%	30%
$A_{\mathrm{eff}}\left[{\mathrm{cm}}^2\right]$	215–1810	20–1200	50–700
FoV $\left[\mathrm{sr}\right]$	2.5	3	3

Equations

$$\frac{d\Phi }{dEd\Omega }=\frac{a⟨\sigma v⟩J}{4\pi m_{\chi }^2}\frac{d{N}_{\gamma }}{d{E}_{\gamma }},$$

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

$$\frac{d{N}_{\gamma }}{dE}=2\delta (E-m_{\chi }).$$

$$\frac{d{N}_{\gamma }}{dE}=\frac{4}{\Delta E}\Theta (E-E_{-})\Theta (E_{+}-E)$$

$$E_{\pm }=\frac{m_{\chi }}{2}\left(1\pm \sqrt{1-\frac{m_{{\pi }^0}^2}{m_{\chi }^2}}\right)$$

$$\frac{d{N}_{\gamma }}{dE}=\delta (E-E_0)+\frac{2}{\Delta E}\Theta (E-E_{-})\Theta (E_{+}-E),$$

$$D(p)=\beta D_0{\left(\frac{p}{p_0}\right)}^{\delta },$$

$$\frac{d\Phi }{dE}=0.013{\left(\frac{E}{1\mathrm{MeV}}\right)}^{-1.8}{e}^{-{\left(\frac{E}{2\mathrm{MeV}}\right)}^2}\mathrm{c}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}\mathrm{s}{\mathrm{r}}^{-1}\mathrm{Me}{\mathrm{V}}^{-1}.$$

$${\varphi }_{\mathrm{bg}}=\sum _{i=1}^4{\theta }_i{\varphi }_i.$$

$${\varphi }_S(E,\Omega |{\vec{\theta }}_S)=\sum _i{\theta }_{S,i}{\varphi }_{S,i}.$$

$${\varphi }_{bg}(E,\Omega |{\vec{\theta }}_{bg})=\sum _i{\theta }_{\mathrm{bg},i}{\varphi }_{\mathrm{bg},i},$$

$${\mathcal{I}}_{ij}(\vec{\theta })=-{⟨\frac{{\partial }^2\ln \mathcal{L}(\vec{\theta }|\mathcal{D})}{\partial {\theta }_i\partial {\theta }_j}⟩}_{\mathcal{D}(\vec{\theta })},$$

$${\varphi }_{\text{bg}}\to (1+\delta (E)){\varphi }_{\text{bg}}.$$

$${\mathcal{I}}_{ij}=\sum _{ab}\frac{{\partial }_i\varphi }{{\varphi }_{\mathrm{bg}}}(E_a)D_{ab}^{-1}\frac{{\partial }_j\varphi }{{\varphi }_{\mathrm{bg}}}(E_b),$$

$$D_{ab}\equiv \frac{{\delta }_{ab}}{\Delta E_a\mathcal{E}(E_a){\varphi }_{\mathrm{bg}}(E_a)}+{\Sigma }_{\text{bg}}(E_a,E_b),$$

$${\Sigma }_{\text{bg}}(E,E^{\prime })=\sum _{k=1,2}{\sigma }_{\mathrm{syst}}^{(k)}(E){\sigma }_{\mathrm{syst}}^{(k)}(E^{\prime }){\rho }_k(E,E^{\prime }),$$

$${\rho }_k(E,E^{\prime })=e^{-\frac{1}{2}{\left(\frac{\ln E/E^{\prime }}{w_k}\right)}^2}.$$

$${\mathcal{I}}_{ij}={\int }_{{\mathrm{E}}_{\min }}^{{\mathrm{E}}_{\max }}dE{\int }_{\mathrm{ROI}}d\Omega \mathcal{E}(E)\frac{{\partial }_i\varphi (E,\Omega ){\partial }_j\varphi (E,\Omega )}{{\varphi }_{\mathrm{bg}}(E,\Omega |\vec{\theta }=\vec{1})}.$$

$${\theta }_{S,i}^{UL}=Z(p)\sqrt{{\mathcal{I}}_{ii}^{-1}},$$

$$p_{\mathrm{ann}}\equiv \frac{f_{\mathrm{eff}}⟨\sigma v⟩}{m_{\chi }}<4.1\times 10^{-31}\mathrm{c}{\mathrm{m}}^3{\mathrm{s}}^{-1}\mathrm{Me}{\mathrm{V}}^{-1}.$$

$$\frac{d\sigma }{d{E}_{\gamma }}=\frac{\pi r_e^2}{m_e{\gamma }^2{\beta }^2}\left(\frac{-\frac{3+\gamma }{1+\gamma }+\frac{3+\gamma }{k}-\frac{1}{k^2}}{{\left(1-\frac{k}{1+\gamma }\right)}^2}-2\right),$$

$$\frac{dN}{d{E}_{\gamma }dt}=N_{e^{-}}\beta c\frac{d\sigma }{d{E}_{\gamma }}.$$

$$\begin{array}{rl}\frac{dN}{d{E}_{\gamma }}& =\frac{1}{{\sigma }_{\mathrm{tot}}}\frac{d{\sigma }_{FSR}}{d{E}_{\gamma }}\\ & =\frac{\alpha }{\pi }\frac{1}{E_{\gamma }}\left[\ln \left(\frac{s^{\prime }}{m_e^2}\right)-1\right]\left[1+{\left(\frac{s^{\prime }}{s}\right)}^2\right],\end{array}$$

$$\begin{array}{rl}\frac{dN}{d{E}_{\gamma }}=\frac{2\alpha }{\pi E_{\gamma }}& \{x^2+2x\left[\mathrm{L}{\mathrm{i}}_2\left(\frac{m_{\varphi }-2m_e}{m_{\varphi }-m_e}\right)-\mathrm{L}{\mathrm{i}}_2(x)\right]+(2-x^2)\ln (1-x)+\\ & \left[\ln \left(\frac{m_{\varphi }^2}{m_e^2}\right)-1\right]\left[2-x^2+2x\ln \left(\frac{\left(m_{\varphi }-m_e\right)x}{m_{\varphi }-2m_e}\right)-\frac{\left(m_{\varphi }^2-2m_e^2\right)x}{\left(m_{\varphi }-m_e\right)\left(m_{\varphi }-2m_e\right)}\right]-\\ & \frac{x}{2m_e^2-3m_{\varphi }{m}_e+m_{\varphi }^2}[2m_e^2\left(2-\ln \left(\frac{m_e^2{x}^2}{(m_{\varphi }-2m_e{)}^2(1-x)}\right)\right)\\ & -3m_e{m}_{\varphi }\left(\frac{4}{3}-\ln \left(\frac{m_e(m_{\varphi }-m_e)x^2}{(m_{\varphi }-2m_e{)}^2(1-x)}\right)\right)\\ & +m_{\varphi }^2\left(1-\ln \left(\frac{(m_{\varphi }-m_e{)}^2{x}^2}{(m_{\varphi }-2m_e{)}^2(1-x)}\right)\right)]\}\end{array}$$

$$D(p)=\left\{\begin{array}[]{lr}D_{0}\left(\frac{p_{b}}{\rho_{0}}\right)^{\delta}\left(\frac{p}{p_{b}}\right)^{\delta_{l}}&\text{ for }p<p_{b}\\ D_{0}\left(\frac{\rho}{\rho_{0}}\right)^{\delta}&\text{ for }p_{b}\leq p,\end{array}\right.$$

Full text

1 Introduction

The past two decades have seen a rapid development of $\gamma$-ray astronomy. The Spectrometer on Integral (SPI) has clearly detected and characterised the diffuse 511 keV electron-positron annihilation signal [1, 2, 3, 4, 5], and the Fermi Large Area Telescope (LAT) has revolutionized the field at energies above $\mathcal{O}(1\mathrm{\,GeV})$ [6]. Yet, sensitivity in the MeV–GeV range has notoriously trailed behind. Existing measurements of diffuse $\gamma$-rays were taken with COMPTEL and EGRET [7, 8]. A number of instruments have been proposed over the last years to cover this so-called ‘MeV gap’ [9]. Currently, the community focuses on various proposals, like e-ASTROGAM111e-ASTROGAM is proposed as ESA M5 mission. [10] and AMEGO222See https://pcos.gsfc.nasa.gov/physpag/probe/AMEGO_probe.pdf. (other previously proposed missions include, e.g., COMPAIR [11] and ADEPT [12]), which could be realized in the mid- and long-term future in the late 2020s. Whatever instrument is ultimately realized, it is expected that it would improve the sensitivity of current dark matter searches with MeV $\gamma$-rays by 2–3 orders of magnitude [13, 14].

Annihilating or decaying dark matter (DM) can leave observable features in the electromagnetic spectrum (for a review see Ref. [15], and references therein). One goal of the proposed MeV telescopes would be to look for signals from particle DM in the diffuse MeV sky. In particular, these telescopes would be sensitive to DM with masses in the MeV–GeV range, henceforth ‘MeV DM’. At this mass scale, the number of kinematically allowed final states becomes limited to neutrinos, pions, photons and light leptons. For DM with masses $\lesssim 10\mathrm{\,MeV}$ in thermal equilibrium in the early universe the strongest constraints come from big-bang nucleosynthesis (BBN), where it can alter the elementary abundances, and from the cosmic-microwave background (CMB) where it affects the effective number of neutrino species. [e.g., 16, 17, 18, 19]. These constraints are independent of whether DM annihilates via $s$- or $p$-wave. In addition, late-time energy injection from DM, between the epochs of recombination and reionization, can leave a strong imprint in the CMB [e.g., 20, 21, 22]. However, since thermal velocities at recombination are low, such bounds are only strong for $s$-wave annihilating DM. Models of DM with MeV masses that can satisfy these constraints and yield the correct relic density include self-interacting DM [23, 24, 25, 26, 27], ‘cannibal’ DM [28] and strongly-interacting DM [29].

Detecting MeV DM is challenging in a number of ways. Direct detection experiments based on nuclear recoils are not sensitive below a GeV, and bounds from ionization are typically weak [30], although some alternative detection principles have the potential to substantially improve current sensitivities in the future [e.g., 13, 31, 32]. Electrons and positrons from MeV DM annihilation or decay are very hard to detect through measurements of the top-of-atmosphere positron fraction because solar modulation cuts-off the cosmic-ray (CR) spectrum at $\mathcal{O}(1\mathrm{\,GeV/nucleon})$ [33]. A notable exception is here Voyager 1, which has left the Heliosphere, and can constrain MeV DM through measurements of the local lepton spectrum [34].

However, indirect detection of the gamma-ray emission from (leptonic) DM annihilation or decay is possible [e.g.  35]. One notable example is the 511 keV anomaly, which has been explained in terms of DM annihilating to $e^{+}e^{-}$-pairs (see e.g. [36] or, for a review, [37] and references therein). However, many models considered in this context (based on light, thermally produced WIMPs) are now ruled out by the early-universe bounds mentioned above  [38]. Thus, a DM explanation of the 511 keV line may require for instance non-thermal or exciting DM [38, 39].

One promising general way to look for MeV DM is to search for spectral features in diffuse $\gamma$-rays. These features can come from the prompt $\gamma$-rays produced in the annihilation, or they are related to secondary emission in leptonic annihilation channels (additional spectral features at sub-GeV energies can come from the decay of meson excited states, as pointed out recently in Ref. [40]). Conservative upper limits on prompt radiation were presented in Refs. [13, 14], without subtraction of any diffuse backgrounds from the data. Moreover, leptonic final states like $e^{+}e^{-}$-pairs can play a particularly interesting role. Aside from the prompt final-state radiation (FSR) signature, there will be a large flux of secondary $\gamma$-rays, dominated by in-flight annihilation (IfA) of positrons. The relevance of IfA was pointed out long ago in Refs. [41, 42, 43], and it has been used to constrain DM explanations of the 511 keV line signal [44, 45]. In addition, there will be a sizeable bremsstrahlung signal. Bremsstrahlung is often ignored in computations of the $\gamma$-ray spectrum from WIMPS, but it becomes increasingly important when going to lower DM masses [46].

In this paper, we study the prospects for a future $\gamma$–ray experiment to detect DM through diffuse $\gamma$–rays originating in the Galactic halo. We will use the characteristics of the proposed e-ASTROGAM [10] as reference. However, we emphasize that our results are also representative for other similar missions, and comment on how the results change with observation time, effective area or energy resolution where appropriate. For the first time, we perform an elaborate study of the detection opportunities of the secondary emission in case of DM annihilation into leptons. In addition, we attempt to derive more optimistic, but also more realistic, projected upper limits for the phenomenologically most interesting final states by modeling the expected uncertainties in the diffuse background. To this end, we assume that remaining systematic uncertainties in these future missions will be of similar size as the ones from the Fermi-LAT today. These uncertainties are then incorporated into a novel statistical approach based on Fisher forecasting [47] (see also Ref. [40]) which we apply to derive projected upper limits. We concentrate here on the a region around the Galactic center, since we know from the Fermi-LAT and H.E.S.S. [48, 49, 50, 51] that this region provides the best (while still reasonably robust) probe for spectral signatures from DM annihilation.

The structure of this paper is as follows: in sections 2 and 3 we discuss the different annihilation channels and their resulting $\gamma$-ray spectra. Sec. 3 is fully devoted to the secondary $\gamma$-rays. We discuss background modeling and the Fisher formalism used to calculate upper limits in Sec. 4. Projected upper limits are presented in Sec. 5 and we end with a discussion and conclusions in Sec. 6 and 7.

2 Dark matter annihilation channels and photon signals

For MeV DM particles, $\chi$, only a few kinematically-allowed two-body annihilation channels exist [14]. We consider here the following processes.

• •

$\chi\chi\to\gamma\gamma$: A photon pair

• •

$\chi\chi\to\gamma\pi^{0}$: A neutral pion and a photon

• •

$\chi\chi\to\pi^{0}\pi^{0}$: Neutral pions

• •

$\chi\chi\to\bar{\ell}\ell$: Light leptons (with $\ell=e,\mu$)

• •

$\chi\chi\to\phi\phi$ and $\phi\to e^{+}e^{-}$: Cascade annihilation

We will here present a complete discussion of the phenomenology of all of these channels in the MeV–GeV regime. For instance, we provide updated prospects for the detectability of DM models in which annihilation proceeds through first-generation quarks into neutral pions or directly into gamma-ray lines as discussed in [14]. However, much of the paper is focused on DM annihilating into charged leptons, since the evaluation the expected signal is more complex than in the other cases. For each channel we will assume for simplicity that the branching ratio is 100%.

2.1 General calculation of gamma-ray signal

The $\gamma$-ray flux resulting from DM annihilation can be split up into two components: a primary and a secondary component.

The primary component is composed of all photons that are produced directly in the annihilation process. Their differential flux is given by [e.g., 52]

[TABLE]

where $\left<\sigma v\right>$ is the velocity-averaged annihilation cross section, $J$ is the astrophysical factor which encapsulates all information about the DM distribution, $m_{\chi}$ is the DM mass, $dN/dE_{\gamma}$ the differential $\gamma$-ray yield per annihilation and $a=\frac{1}{2}\,(\frac{1}{4})$ if DM is (is not) self-conjugate. Throughout this work we assume that DM is self-conjugate.

The so-called astrophysical-, or J-factor, is given by

[TABLE]

where $\rho(r)$ is the DM density as function of the Galacto-centric radius $r$. We adopt here a Navarro-Frenk-White (NFW) density profile [53], $\rho_{\mathrm{NFW}}(r)=\rho_{0}/((r/r_{s})^{\gamma}(1+r/r_{s})^{(3-\gamma)})$, with a local DM density of $\rho_{\odot}=0.4\mathrm{\,GeV\,cm^{-3}}$ [e.g.  54, 55, 56, 57, 58, 59], scale radius $r_{s}=20\mathrm{\,kpc}$ and slope $\gamma=1$. Since DM annihilation scales as $\rho^{2}$, our results are – as typical for searches in the inner Galaxy – very sensitive to the adopted profile. We will address how different assumptions on the profile affect our main results in Sec. 6.

Regarding the DM annihilation cross section, given the usual decomposition $\left(\sigma v\right)=a+bv^{2}$, with $a$ the $s$-wave and $bv^{2}$ the $p$-wave contribution, we require that the $p$-wave term sets the relic density and still dominates at the time of recombination, in order to avoid the CMB constraint from late-time energy injection mentioned in the introduction. For this reason, we expect a value much below the canonical $3\cdot 10^{-26}\mathrm{\,cm^{3}\,s^{-1}}$. For example, for a thermal relic purely annihilating through $p$-wave processes, the expected annihilation cross section times velocity in the Galaxy today is of the order $\sim 10^{-31}\mathrm{\,cm^{3}\,s^{-1}}$ (see discussion in Sec. 5.3).

The prediction of the secondary $\gamma$-ray flux from charged particles (here electrons) is significantly more complex. The $\gamma$-ray spectrum will depend on the energy-losses of the electrons, and on the radiative process underlying the emission. In addition, the morphology will no longer trace the DM squared distribution directly, but rather depend on how far and in what direction the electrons have propagated. These environmental impacts are discussed extensively in section 3.

2.2 Gamma-ray lines and neutral pions

We briefly summarize the analytical expressions for the various prompt annihilation spectra that we consider in this work.

1. (i)

$\chi\chi\rightarrow\gamma\gamma$. For this channel, the photon spectrum per annihilation is given by

[TABLE] 2. (ii)

$\chi\chi\rightarrow\pi^{0}\pi^{0}$. The neutral-pion channel leads to a box-shaped $\gamma$-ray spectrum [60, 14],

[TABLE]

per annihilation. Here, $\Theta$ is the Heaviside step function,

[TABLE]

are the kinematic edges of the box, and $\Delta E\equiv E_{+}-E_{-}=\sqrt{m_{\chi}^{2}-m_{\pi_{0}}^{2}}$ denotes the box width. 3. (iii)

$\chi\chi\rightarrow\pi^{0}\gamma$. This channel leads to a box and $\gamma$-ray line, with slightly different kinematics from the discussion above [14].

The prompt photon spectrum per annihilation is

[TABLE]

where $E_{0}=m_{\chi}-\frac{m_{\pi_{0}}^{2}}{4m_{\chi}}$, $\Delta E=m_{\chi}-\frac{m_{\pi_{0}}^{2}}{4m_{\chi}}$, and $E_{\pm}=\frac{m_{\chi}}{2}\left[\left(1+\frac{m_{\pi^{0}}^{2}}{4m_{\chi}^{2}}\right)\pm\left(1-\frac{m_{\pi^{0}}^{2}}{4m_{\chi}^{2}}\right)\right]$.

Finally, we briefly comment on the charged-pion channel, that we neglect otherwise in this work. The dominant decay channel for charged pions is $\pi^{+}\rightarrow\mu^{+}\nu_{\mu}$ [61]. The muon will subsequently decay into an electron and two neutrinos. Therefore, the final electron spectrum will receive two boosts, one from the annihilation and one from the subsequent decay. In addition, FSR can be produced in both the annihilation as well as in the decays. Because this channel does not add any new spectral signatures, but rather smears out the ones studied in the other channels, we ignore the charged pion channel in the present analysis.

2.3 Leptonic annihilation

In the case of annihilation into leptons, we consider three benchmark channels for the DM annihilation channels. We comment here briefly on the role of prompt photons for these channels and leave a detailed discussion about secondary emission for the next section.

1. (i)

$\chi\chi\rightarrow e^{+}e^{-}$. In this case electrons and positrons are injected mono-energetically. There is a significant contribution from final-state radiation (FSR), which we model following [62]. An analytic expression is provided in appendix A. 2. (ii)

$\chi\chi\rightarrow\phi\phi$, $\phi\rightarrow e^{+}e^{-}$. Here, DM annihilates via a (scalar) mediator which subsequently decays into an electron/positron pair. Now the positron spectrum will be boosted and has a box-like shape333Under the assumption that $m_{\mathrm{e}}\ll m_{\phi}<m_{\chi}$ and $\phi$ is a scalar. [e.g.  63, 60]. An analytic expression for FSR in such models is given in Eq. (6) of [64] and for completeness repeated in appendix A.

FSR for this class of models is suppressed, enhancing the relative importance of secondary emission. On the other hand, the final electron spectrum is boosted, softening any spectral features. 3. (iii)

$\chi\chi\rightarrow\mu^{+}\mu^{-}$. Like in the cascade channel, the injected electron spectrum gets boosted. On the other hand, this channel is accompanied by a large amount of FSR, since FSR can be produced at two stages: when dark matter annihilates to muons and in the subsequent decay of the muon to an electron and two neutrinos. Both the positron and FSR spectra are obtained from DarkSusy [65]444P. Gondolo, J. Edsjö, P. Ullio, L. Bergström, M. Schelke, E.A. Baltz, T. Bringmann and G. Duda, http://www.darksusy.org.

Finally, we note that we do not consider the channel $\chi\chi\rightarrow\tau^{+}\tau^{-}$, since the tau lepton is heavy, $m_{\tau}=1.78\mathrm{\,GeV}$, and therefore only relevant for a very small window in our considered energy range.

3 Gamma-ray signals from secondaries

3.1 Radiative processes, cooling and timescales

As mentioned above, we concentrate in this analysis on the annihilation signal from the inner Galaxy (projected limits from dwarf spheroidal galaxies were discussed in Ref. [14]). The higher statistics from this region makes it easier to identify spectral features, which are the main focus of our work. Furthermore, secondary emission components are in most cases expected to be stronger. However, in order to correctly predict $\gamma$-rays from secondaries it is essential to understand all cooling and radiative processes. We begin this section with a brief discussion of all relevant processes and timescales.

The relevant radiative processes in our energy window are bremsstrahlung, inverse-Compton scattering (ICS) and synchrotron radiation. The first two yield X-ray and $\gamma$-ray photons, whereas the latter leads to radio emission. In addition, ionization and Coulomb losses are important contributors to the overall energy losses. Emissivities and energy-loss timescales are well established and can be found in Refs. [66, 33, 67].

In addition, we model in-flight annihilation (IfA) of positrons. IfA results from positrons colliding with ambient electrons. The cross section for this process is implemented following Refs. [68, 44]. An expression for the differential cross section and resulting photon flux for a single positron is given in appendix A. Note that at very low energies, below those considered in this work, positrons can also annihilate through the formation of positronium, a bound state of a positron and an electron. Positronium annihilation results in monochromatic photon emission at 511 keV from the singlet (p-Ps) state and continuum emission from the triplet (o-Ps) state.

In Fig. 1 we show the energy-loss and annihilation timescales for electrons and positrons from 1 MeV to 1 TeV in a medium comparable to that of the Galactic bulge. As can be seen, at energies $\lesssim 100\mathrm{\,MeV}$ IfA starts to dominate over ICS and synchrotron losses. Also bremsstrahlung becomes more important. Energy losses are dominated by ionization and Coulomb losses, with their relative importance depending on the ionization fraction of the medium. Note that bremsstrahlung, IfA, ionization- and Coulomb losses all depend on the gas densities, whereas ICS depends on the interstellar-radiation field and synchrotron emission on the magnetic field.

3.2 Cosmic-ray propagation

We model the transport of CR electrons and positron in the Galaxy using the numerical code DRAGON [69].555www.dragonproject.org. DRAGON is designed to simulate all the relevant processes related to Galactic CR propagation, in particular: diffusion, reacceleration, convection, cooling (due to synchrotron, bremsstrahlung, inverse-Compton, Coulomb and ionization), catastrophic losses (annihilation), and spallation.

The transport equation is solved for all CR species, from heavy nuclei down to protons, antiprotons, and leptons. The code implements both the nuclear spallation network taken from the public version of GALPROP666GALPROP project web page, http://galprop.stanford.edu (see e.g. [70] and references therein), and a complete cross-section database obtained with the FLUKA code [71]. The code works in 2D and 3D mode; in the following we assume azimuthal symmetry and work in the two-dimensional mode.

The annihilation and energy-loss rate is computed adopting the following astrophysical ingredients:

• •

The gas distribution is based on an azimuthally-averaged implementation of the detailed three-dimensional model described in Ref. [72] for the Galactic bulge ($R<3$ kpc), while the model presented in Ref. [73] is used for the larger Galactocentric radii. Within $R\lesssim 0.7$ kpc the ionization fraction is $n_{\mathrm{H^{+}}}/(n_{\mathrm{H}}+2n_{\mathrm{H_{2}}})\approx 0.1$.

• •

The magnetic field is taken from [74] and is derived by a wide set of Faraday Rotation measurements.

• •

The interstellar radiation field is taken from the public version of GALPROP; it is described in [75, 76].

Diffusion of sub-GeV electrons remains largely unconstrained (partially due to the lack of observational data). We assume here isotropic diffusion with power-law dependence on momentum, $p$. The diffusion coefficient is

[TABLE]

where $\beta=v/c$, $p_{0}=4\mathrm{\,GeV}$, $\delta=\frac{1}{2}$ and $D_{0}=4\times 10^{28}\mathrm{\,cm^{2}\,s^{-1}}$. We verified with DRAGON that these parameters are compatible with the current AMS-02 measurements of the boron-over-carbon ratio in the GeV-TeV range [77]. We comment on the effects of changing these assumptions on the diffusion coefficient in Sec. 6. For simplicity, we do not consider diffusive reacceleration in this work.

Soft X-rays provide evidence for a galactic wind [78]. Following Ref. [78] we implement a cylindrical Galactic wind of radius $R_{\mathrm{wind}}=3\mathrm{\,kpc}$ and velocity $v_{\mathrm{wind}}=250\mathrm{\,km\,s^{-1}}$ centered on the Galactic center into our benchmark model. As an alternative we also show results in absence of a wind in appendix B.

The relevant timescales for diffusion and advection are included in Fig. 1. Under the above benchmark assumptions, advection dominates over diffusion below $\sim 10\mathrm{\,GeV}$. Moreover, cooling is much faster than diffusion. In the absence of a wind, sub-GeV particles would approximately cool in-situ. However, when a wind is present ($v_{\mathrm{wind}}=250\mathrm{\,km\,s^{-1}}$), a particle at the Galactic Center (GC) will be advected over $5^{\circ}$ on the sky. We apply a region-of-interest of $10^{\circ}\times 10^{\circ}$ centered at the GC, implying that in the presence of a wind a particle is more likely to be advected out of our ROI than to radiate its energy.

3.3 Secondary gamma-rays

Next, the $\gamma$-ray spectrum is obtained by convolving the steady-state electron and positron spectra with either the gas-distribution (bremsstrahlung and IfA) or the interstellar-radiation field (ICS) and performing an integral over the line-of-sight. We ignore emission from synchrotron radiation, since it will be in radio, but we do include it as a cooling term. In Fig. 2 we show an example $\gamma$-ray spectrum for a $m_{\chi}=60\mathrm{\,MeV}$ DM particle annihilating to $\chi\chi\rightarrow e^{+}e^{-}$. Note that there is a kinematic cutoff in the IfA spectrum at $E_{\gamma}=m_{e}/2$ [79] (for more details see appendix A).

4 Background modeling and sensitivity projections

4.1 Instrumental details

We assume a future instrument similar to the mission concept e-ASTROGAM [10]. Marginal differences put aside, our results also roughly apply to other proposed mission such as AMEGO, COMPAIR [11] or ADEPT [12] within their energy range. An important reason for this is that our projected sensitivity turns out the be systematics limited over most of the considered DM mass range. This is a consequence of our use of a relatively large ROI on the inner Galaxy. This makes slight variations in the effective area irrelevant for the projected limits. On the other hand, a better spectral resolution, such as that of COMPAIR, could potentially improve the projected limits in the pair–regime by a factor two (see Tab. 1). We provide an overview of instrumental details and compare instruments in Tab. 1. These instruments can fill the MeV gap and as we will show have great prospects to hunt for MeV DM. At low energies ($<10\mathrm{\,MeV}$) the detection principle is based on Compton scattering (similar to COMPTEL) and at higher energies on the detection of $\gamma$-rays via pair-conversion (like Fermi-LAT).

All instrumental details applied in our analysis are taken from Fig. 17 and tables 2 and 3 of Ref. [10]. The effective area ($A_{\mathrm{eff}}(E)$) is modeled as a function of energy and we interpolate between the reference values. It is typically of the order $A_{\mathrm{eff}}\sim\mathrm{few}\times 100(1000)\rm\,cm^{2}$ in the Compton (pair-production) domain. The spectral resolution is modeled as a Gaussian with energy dependent width. It is of the order $\Delta E/E\sim 0.01(0.3)$, in the Compton (pair-production) domain. The excellent energy resolution in the Compton domain yields great power to search for sharp spectral features such as $\gamma$-ray lines, as we will demonstrate below (Fig. 3). Finally, we adopt an effective observation time on our ROI of $T_{\mathrm{obs}}=1\mathrm{\,yr}$. Considering that the field-of-view of e-ASTROGAM is $2.5\mathrm{\,sr}$, this corresponds to $\sim 5$ years of pure observation time with uniform sky coverage.

4.2 Background components

Our diffuse background model consists of four components as shown by the black lines in Fig. 2. Three of them – the bremsstrahlung, $\pi^{0}$ and inverse-Compton components – characterize the diffuse background at higher energies. We take the spectral templates from [80] (see also [81]). Those models were computed with the GALPROP code [82], and fitted to data in the window $|l|<30^{\circ},|b|<10^{\circ}$. Since the region of interest (ROI) for this analysis is $|l|<5^{\circ},|b|<5^{\circ}$, we rescale those models to match our ROI, assuming the same morphology computed with DRAGON at 1 GeV. We checked that the overall $\gamma$-ray intensity of our background model in the ROI is compatible with Fermi-LAT data.

The fourth component characterizes the background at lower energies and is modeled as a power-law, following [44], and including a super-exponential cutoff:

[TABLE]

This component, possibly originating from Inverse-Compton scattering, is introduced to get a reasonable agreement with the COMPTEL data in the region $|l|<30^{\circ},|b|<5^{\circ}$ [44].

We will henceforth refer to this emission as $\mathrm{ICS_{lo}}$ and to the inverse-Compton model at higher energies as $\mathrm{ICS_{hi}}$. Similarly to the other three background models, we map the intensity to that in our ROI utilizing the DRAGON ICS template.

Our full background model is then described by

[TABLE]

Here $\phi(E,\Omega)=d\Phi/dEd\Omega$ and $i=1,2,3,4$ refers to bremsstrahlung, $\mathrm{ICS_{hi}}$, $\pi^{0}$ and $\mathrm{ICS_{lo}}$, respectively. The parameters $\theta_{i}$ are the normalization of the various components. The baseline background model has $\theta_{i}=1$ for components $i=1,2,3$ and for component $i=4$ such that the latter reproduces Eq. 4.1 within $|l|<30^{\circ},|b|<5^{\circ}$

It is important to emphasize that – for the purposes of this paper – an approximate treatment of the background suffices. A dedicated, more accurate modeling of the diffuse $\gamma$-ray sky from MeV to GeV and a corresponding, consistent modeling of the dark matter signal is well beyond the scope of the present paper and will be left for future work.

4.3 Projected limits from Fisher forecasting with correlated background

systematics

We compute here projected 95% confidence level (CL) upper limits on any DM signal, assuming that no signal is present in the data. Fisher forecasting is applied to calculate projected limits for a mission similar to e-ASTROGAM [47].

We consider a signal spectrum given by the differential flux $\phi(E,\Omega)=d\Phi/dE(\Omega)$ in units $\mathrm{ph\,cm^{-2}\,s^{-1}\,GeV^{-1}}$, composed of various additive components with subscript $i$,

[TABLE]

In our case, we include the components $i\in\{\mathrm{bremsstrahlung,\,FSR,\,ICS,\,IfA}\}$.

Additionally, we consider the background and foreground components

[TABLE]

as described in Eq. 4.2. The total emission is then given by $\phi=\phi_{S}+\phi_{\mathrm{bg}}$.

The Fisher information matrix is an $N\times N$ matrix where $N$ is the number of parameters $\vec{\theta}$. Given a likelihood function, $\mathcal{L}(\vec{\theta}|\mathcal{D})$, it is defined as

[TABLE]

where the average is taken over multiple realizations of the data, $\mathcal{D}$. Practically, the average is taken over $\mathcal{D}(\vec{\theta}_{\mathrm{S}}=\vec{0},\vec{\theta}_{\mathrm{bg}}=\vec{1})$ when setting upper limits.

From here on we will assume an unbinned Poisson likelihood for the description of mock data. Additionally, we include an effective model for correlated instrumental systematics and/or background model uncertainties in our analysis, by making the substitution

[TABLE]

Here, $\delta(E)$ parametrizes general fractional deviations from the nominal background model, which are assumed to be correlated over smaller and larger energy ranges, as specified below. One can then show that the Fisher information becomes [47]

[TABLE]

with $\partial_{i}\phi\equiv d\phi/d\theta_{i}$. Here, $E_{a}$ (and equivalently $E_{b}$) refer to a dense grid of reference energies at which $\delta(E)$ and all other fluxes are evaluated. It has to be sufficiently dense to capture all spectral variations. For this analysis we use 2000 logarithmically spaced bins between $0.5\mathrm{\,MeV}$ and $5\mathrm{\,GeV}$. We note that the secondary $\gamma$-ray  spectra obtained through DRAGON have a more coarse binning due to computational constraints. In practice we interpolate over these spectra to get a sufficiently fine binning.

Furthermore, $D_{ab}$ is defined as

[TABLE]

where $\Delta E_{a}$ is the energy step between two consecutive values of $E_{a}$. Furthermore, $\mathcal{E}(E)\equiv T_{\mathrm{obs}}A_{\mathrm{eff}}(E)$ denotes the exposure,777In principle, the exposure is a function of not only energy, but also sky-position through $T_{\mathrm{obs}}(\Omega)$. However, for this paper we assume uniform sky coverage. which is an energy dependent quantity. Finally, $\Sigma_{\text{bg}}$ denotes the covariance matrix that describes the background model uncertainties encoded in $\delta(E)$ (not to be confused with the covariance matrix of the model parameters, which equals $\mathcal{I}^{-1}$). It can be thought of as the covariance matrix of a Gaussian random field, with $\Sigma_{\text{bg}}(E,E^{\prime})=\langle\delta(E)\delta(E^{\prime})\rangle$, where the average is taken over many realizations of the function $\delta(E)$. It is here parameterized as

[TABLE]

where we assumed two independent contributions. Here $\sigma^{(k)}_{\mathrm{syst}}(E)$ is the overall magnitude of the systematic uncertainty at energy $E$, while $\rho_{k}(E,E^{\prime})$ parametrizes the correlation between systematics at different energies and equals one along the diagonal. We adopt the simple form (motivated by a log-normal distribution)

[TABLE]

As mentioned above, we use here two independent covariance matrices, which are added ($k=1,2$), and which correspond to background systematic with short and long correlation scales in energy space. We adopt 2% ($\sigma^{(1)}_{\text{syst}}=0.02$) as the short-scale correlated systematic, with a correlation length of just $1\%$ ($w_{1}=0.01$) in energy space, which is motivated by the results from Refs. [83, 50] and presumably mostly of (unknown) instrumental origin. For systematics correlated over a larger energy range, possibly related to diffuse emission modeling and uncertainties in the effective area, we assume a value of 15% ($\sigma_{\text{syst}}^{(2)}=0.15$), with the correlation length taken as 0.5 dex ($w_{2}=1.15$). This is reasonably representative of what is found when analyzing Galactic diffuse emission along the Galactic disk [84].

We note that our treatment of systematic uncertainties for future instruments remains necessarily uncertain, although we believe that our choices are plausible. However, we checked that our qualitative results are not sensitive to changes in in these parameters, and the quantitative results are fairly independent of $w$ and behave as expected for changes in $\sigma$. For a more elaborate discussion see Sec. 6.

When systematic uncertainties are absent, $\Sigma_{\text{bg}}\rightarrow 0$, the Fisher information can be written in the more common form in terms of the signal-to-noise ratio [47],

[TABLE]

Finally, in the background-limited regime we can construct an upper limit of $100\cdot(1-p)\%$ confidence level (CL) on a signal parameter, $\theta^{UL}_{S,i}$, by inverting the Fisher-information,

[TABLE]

where $Z(p)=\mathrm{ppf}(1-p)$, with $\mathrm{ppf}(x)$ the inverse of the cumulative-distribution functions, also known as the percent-point function. For 95% CL (or $p=0.05$) $Z(p)=1.645$.

5 Results

5.1 Projected limits for $\gamma$-ray lines and pions

Figure 3 shows the projected upper limits from the Galactic center for annihilation into mono-chromatic photons (red), a neutral pion and a photon (green) and two neutral pions (magenta). All projected limits are 95% CL, for one year of effective exposure and for instrumental specifications similar to those of e-ASTROGAM. The emission is prompt and therefore traces exactly the DM distribution.

The limits from monochromatic photons are only slightly stronger than those from box-spectra due to pion decay. This is a consequence of the spectral resolution in the pair-production domain of e-ASTROGAM, assumed to be $\sim 30\%$, which significantly broadens the line feature: as a result, the line is not much sharper than the box-like feature. The fact that pion decay yields two photons, leading to twice as many photon in $\chi\chi\rightarrow\pi^{0}\pi^{0}$ compared to direct annihilation into photons, further reduces the difference. Below 10 MeV the limits for monochromatic photons improve by an order of magnitude, because of the better spectral resolution in the Compton domain (a factor $\sim 10$ better with respect to the pair–creation regime).

A MeV mission similar to the proposed e-ASTROGAM could outperform CMB constraints for $s$-wave annihilating DM (shown for $\chi\chi\rightarrow\gamma\gamma$ in Fig. 3) by more than one order of magnitude below 1 GeV.

Our projections are compared to Fermi–LAT limits from the inner-Galaxy on monochromatic photons resulting from DM annihilation with $m_{\chi}>200\mathrm{\,MeV}$ as obtained by [50]. The limits portrayed correspond to the analysis optimized for a NFW density profile.

In addition, we compare our limits to the existing and forecasted limits from Ref. [14]. The existing limits, shown as the light-grey shaded regions, are derived from the diffuse $\gamma$-ray flux measured by COMPTEL ($|b|>30^{\circ}$), EGRET ($20^{\circ}<|b|<60^{\circ}$) and Fermi–LAT ($|b|>20^{\circ}$). Our projected limits for the diffuse $\gamma$-ray sky suggest that these constraints can be improved by at least two orders of magnitudes. We note that a reanalysis of existing data and modeling of background could improve the current situation already now significantly, and would be an interesting endeavour.

For the forecast, Ref. [14] assumes an ADEPT–like instrument with a spectral resolution of $\Delta E/E=15\%$ (a factor two better than what we assume), an effective area of $A_{\mathrm{eff}}=600\mathrm{\,cm^{2}}$ and a systematic uncertainty in the backgrounds of 15%. The total observation time of the experiment is set to 5 years. Ref. [14] considers two targets, and the corresponding projections are shown in Fig. 3: the optimistic projection for diffuse $\gamma$-rays above $|b|>30^{\circ}$ (dashed) and for the Draco dwarf-spheroidal galaxy (dotted). A detailed discussion of the differences will be done in Sec. 6.

5.2 Projected limits for leptonic channels

In Fig. 4 we show the spectral constraints (95% CL) obtained for our three reference leptonic channels. We adopt the transport model described in section 3.2: The dominant effect in the low-energy range considered in this work is advection, caused by a Galactic wind modeled with $v_{\mathrm{wind}}=250\mathrm{\,km\,s^{-1}}$ (see also Fig. 1). The impact of diffusion is negligible at low energies for our benchmark scenario, in which we consider a power-law extrapolation of the diffusion coefficient tuned on GeV boron-over-carbon AMS data.

The three panels in Fig. 5 show the ratio of the limit obtained from a single emission component (i.e., bremsstrahlung, FSR, ICS or IfA) over that of the full DM spectrum for the electron/positron, cascade and muon channel, respectively. The plots show the importance of each emission component. In case of mono-energetic injection of electron-positron pairs, in-flight annihilation of positrons dominates the bounds below $\sim 20$ MeV. For the cascade annihilation scenario the upper limits arise predominantly from IfA below $\sim 50$ MeV, since FSR is suppressed. However, the overall limit is somewhat weaker due to the softening of the injected lepton spectrum (see Fig. 8(b) in appendix C). From $\sim 50$–$200$ MeV bremsstrahlung provides the dominant signal. The muonic channel is most easily detectable through FSR at all DM masses. In this case FSR arises at two stages, when DM annihilates to muons, and in the subsequent decay of the muon. ICS dominates the bounds above a few hundred MeV for the direct channel into $e^{+}e^{-}$ and for the cascade channel. However, unlike FSR, IfA and to some extent bremsstrahlung the ICS spectrum is not very peaked and will therefore be more difficult to distinguish from any astrophysical background.

In light-grey we show the recent limits from Voyager on $\chi\chi\rightarrow e^{+}e^{-}$ (Fig. 4 [34]). Future diffuse $\gamma$-ray studies can surpass these limits for MeV DM and cover a broader mass region in general. In addition, we show existing diffuse $\gamma$-ray limits [13] (based on INTEGRAL and COMPTEL data, and without performing any background subtraction). A future dedicated MeV $\gamma$-ray experiment can improve these limits by two to three orders-of-magnitude (less impressive but still interesting improvements could be obtained by a dedicated reanalysis of available data, including a proper modeling of backgrounds).

Figure  4 shows that the CMB limits from Planck (blue solid line) for $\chi\chi\rightarrow e^{+}e^{-}$ through an $s$-wave process are stronger than what is attainable with a MeV mission like e-ASTROGAM. However, as mentioned above thermal relic models will require $s$-wave suppression in order to evade the CMB bounds. In this case the CMB constraints disappear, whereas the limits from diffuse $\gamma$-rays remain.

5.3 Other constraints

We here briefly review limits from the CMB, BBN and the $p$–wave thermal relic cross section to which we compare our results for the projected upper limits from future diffuse $\gamma$-ray studies.

Regarding the CMB bounds for late time energy injection from Planck, these apply to $s$–wave annihilation. The constraint on the so-called annihilation parameter, $p_{\mathrm{ann}}$ is (for the Planck TT, TE, EE+lowP data) [21]:

[TABLE]

Here $f_{\mathrm{eff}}$ is an efficiency parameter that relates the total injected energy to the energy that is used to increase the ionization fraction. We adopt values for $f_{\mathrm{eff}}$ for $\chi\chi\rightarrow\gamma\gamma$ and $\chi\chi\rightarrow e^{+}e^{-}$ from Ref. [22]888Tables available from http://nebel.rc.fas.harvard.edu/epsilon.. Constraints for the cascade channel, which are not shown, are only marginally different from those for the direct annihilation channel [85]. For the muonic channel they are weaker by a factor of a few [22, 85]. For a future experiment that is cosmic-variance (CV) limited the bound can be as strong as $p_{\mathrm{ann}}<8.9\times 10^{-32}\mathrm{\,cm^{3}\,s^{-1}\,MeV^{-1}}$, about five times stronger than the current Planck bound[20]. CMB limits from Planck for late time energy injection are displayed as a blue solid line in Figs 3 and 4. We note that any freeze-out thermal relic DM candidate will require velocity suppressed annihilation in order to evade the CMB bounds. For $p$-wave annihilation the CMB bounds deteriorate by a few orders of magnitude [86].

In addition, for low mass thermal dark matter, there exist constraints from big-bang nucleosynthesis (BBN) and from the effective number of neutrinos as inferred from the CMB by Planck. The latter provide the stronger constraints, ruling out a thermal dark matter candidate in the form of a Majorana (Dirac) fermion below $m_{\chi}<3.5\,(7.3)\mathrm{\,MeV}$, independent of whether the annihilation is $s$- or $p$-wave [17].

Finally, we comment on how the projected upper-limits compare to a typical annihilation cross section for a pure $p$–wave annihilating thermal relic in the Galactic center. For this estimate we assume that the DM velocity distribution can be described by a Maxwell-Boltzmann distribution in the early universe. At freeze-out ($m_{\chi}/T\approx 20$) the DM velocity dispersion is $\sigma_{v,\mathrm{\,fo}}\approx 0.4c$, where $c$ is the speed of light. Today, the velocity dispersion in the inner Galaxy is $\sigma_{v}\approx 4\times 10^{-4}c$ [87]. So $\left<\sigma v\right>=\left<\sigma v\right>_{\mathrm{fo}}\left(\sigma_{v}/\sigma_{v,\mathrm{\,fo}}\right)^{2}\approx 6\times 10^{-32}\,\mathrm{\,cm^{3}\,s^{-1}}$. Here $\left<\sigma v\right>_{\mathrm{fo}}=6\times 10^{-26}\mathrm{\,cm^{3}\,s^{-1}}$ is the $p$-wave annihilation cross section at the time of freeze out [86]. If DM annihilates to monochromatic photons, this cross-section can be probed for $m_{\chi}\lesssim 10$ MeV. On the other hand, the projected limits for leptonic final states are at least one order of magnitude above this reference thermal $p$–wave cross-section.

However, we note that in the vicinity of the supermassive black hole Sagittarius $\mathrm{A^{*}}$ (Sgr $\mathrm{A^{*}}$) where much of our signal arises, the velocity dispersion could potentially be an order of magnitude larger, increasing the annihilation cross section around the black hole by two orders of magnitude [see e.g.  88, 89], potentially bringing a thermal relic that annihilates through pure $p$–wave processes within reach.

6 Discussion

In the previous section we have shown that – under reasonable assumptions for the systematic uncertainties for the backgrounds/foreground ($\sim 2\%$–$15\%$) and for a reasonable benchmark cosmic-ray transport model – future MeV $\gamma$-ray missions have the capability to either detect a signature of thermal relic DM in the case of $p$-wave annihilation (monochromatic photons), or set stringent bounds that are competitive with current CMB bounds for $s$-wave annihilation, and much stronger in the case of $p$-wave annihilation. The existing current indirect detection limits on MeV dark matter can be improved by two to three orders of magnitude in case of diffuse $\gamma$-rays, and by a factor of a few in the case of direct measurements of positrons.

Let us now compare the projected bounds derived in this work with previous studies, comment on the impact of our assumptions and discuss a potential background from nuclear de-excitations.

Comparison to previous works. As mentioned in Sec.5, our projections imply that future experiments can improve on existing limits from Refs. [14, 13] by 2–3 orders of magnitude. This difference is partially due to the fact that we consider a better region of interest (the inner Galaxy), and a more powerful instrument, but mostly because existing limits were derived without any background subtraction. In the leptonic annihilation case, another difference between this work and Ref. [13] is that we consider all signal components, both primary and secondary. This is very important since, in the the inner Galaxy, secondaries can actually provide the dominant $\gamma$-ray  signal. However, it should be pointed out that Ref. [13] considers higher Galactic latitudes, where in general gas densities are lower and therefore timescales for gas-related radiative processes longer. As such, a wind or enhanced diffusion can more easily reduce the secondary signal and enhance the importance of FSR.

When comparing our projections, based on inner Galaxy observations, to those from Ref. [14], we find that our results only differ by a factor of a few from their projected limits for the Draco dwarf-spheroidal galaxy. This is merely coincidental. The ratio of integrated J-factors between Draco and our ROI is $J_{\mathrm{inner\text{--}Galaxy}}/J_{\mathrm{Draco}}\sim 10^{5}$. However, the expected background in a dwarf galaxy is only a handful of photons, allowing for very strong limits despite a limited signal. Their ‘diffuse’ projections are however two orders of magnitude weaker. We believe this difference can be attributed to significant differences in the ROI as well as different assumptions on the expected instrumental and background systematics. More specifically, Ref. [14] assumes a 15% systematic uncertainty on the backgrounds per bin, which is similar to our value for uncertainties with a large-range correlation in energy space. However, we adopt only a 2% systematic uncertainty for fluctuations that have the width of $\gamma$–ray lines, as motivated by previous Fermi-LAT searches [83, 50].

Assumptions on systematics. As outlined in Sec. 4, we include systematics in our Fisher forecast and assume a 2% uncertainty on small-scales (1% in energy), and 15% on large scales (corresponding to 0.5 dex in energy). These parameters are chosen such that they mimic our understanding of systematics with Fermi–LAT today. We checked that varying the correlation lengths of the small- and large-scale correlations only has a small impact on our results. As mentioned above, the projected sensitivity to DM annihilation signals in this work are systematics limited. Consequently, doubling the observation time only yields minor ($\sim 10$–$20\%$) improvements in the limits.

Cosmic-ray propagation. Our secondary signal depends on the propagation model of MeV leptons. Understanding CR propagation in this low-energy domain is extremely challenging, and little or no constraints arise from available data. In the simplest benchmark case, we model the diffusion coefficient as a single power law in rigidity and extrapolate down from the GeV domain where secondary-to-primary ratios set the normalization and slope. In this case, diffusion hardly impacts MeV DM (see Fig. 1), and CR transport is dominated by the competition between advection and energy losses.

We also consider a different scenario characterized by a significantly faster diffusion below 1 GeV. A viable physical picture behind this scenario is e.g. the dissipation of the magneto-hydrodynamic waves responsible for CR confinement due to the resonant interaction with low-energy CRs themselves (see e.g. [90]). However, even if the diffusion coefficient flattens at rigidity of $\sim 1$ GeV due to this effect, our results are not severely affected, since cooling and advection timescales are still faster in the energy range we are interested in.

In the absence of a wind, the secondary spectra for low DM masses will be enhanced by a factor of a few, since particles no longer escape their injection site, but instead cool and radiate in-situ. The projected limits for this scenario increase marginally, since secondary emission becomes yet more important. For a comparison of the projected limits with and without wind see Appendix B. Note that our benchmark model that includes an advective wind of $250\mathrm{\,km\,s^{-1}}$ starting from the plane appears conservative, since the wind velocity, in particular close to the plane, could be much smaller [91], thereby increasing the advection timescale. However, as argued above, our overall results do not critically depend on the exact transport parameters chosen.

$\gamma$–ray lines from nuclear de-excitations. A potential background that can produce a sharp feature in the low-energy $\gamma$-ray spectrum and mimic a DM annihilation signal can come from the inelastic interactions of low-energy cosmic rays with interstellar gas. This process can produce a rich spectrum of de-excitation nuclear $\gamma$-ray  lines from $\sim 0.1$ to $\sim 10$ MeV, spatially correlated with the gas distribution [92]. For this reason, a careful morphological study of a possible future signal will be needed in order to disentangle astrophysical and DM interpretations.

DM profile. We briefly comment on our choice of dark-matter profile. This work applies a standard NFW profile with an inner slope of $\gamma=1$. However, the actual density profile of the Milky–Way is still uncertain as its determination requires detailed knowledge of the baryonic-mass distribution [e.g.  93, 94]. Since annihilation scales as the density squared, having a different profile, such as the Einasto profile [95] or a cored profile, can lead to relevant differences. For the prompt signal, the limits will simply scale with the difference in J-factor in our ROI, which for Einasto is $1.5\times$ larger and for a cored-isothermal profile ($r_{\mathrm{core}}=4.4\mathrm{\,kpc}$) $10\times$ smaller [52]. For the secondaries this dependence is not one-to-one, due to the effects of diffusion, cooling and radiation, of which the latter two depend on the ambient medium. However, no large deviations are expected, in particular since particles are quickly moving away from the densest gas region in our analysis due to the wind.

511 keV emission. Finally, let us briefly comment on 511 keV emission. When positrons are injected, there will be 511 keV line emission. However, for the projected upper limits on the cross section, the associated 511 keV signal is only a fraction of the total line emission. This should be obvious, since it was already pointed out by Ref. [44] that in order to explain the full 511 keV signal with DM there should be an observable in-flight annihilation signal. Since there are many candidates to produce the low-energy positrons that produce the 511 keV signal [see 37, for a review], the uncertainty in what fraction could be due to DM is large. Under the assumption that most 511 keV emission is astrophysical in origin, it is virtually impossible to detect DM through this channel in the Inner Galaxy. On the other hand, we have demonstrated that both prompt and secondary emission resulting from DM annihilation can produce spectral signatures that should be detectable in diffuse $\gamma$-rays where the fore- and backgrounds are better understood. Therefore, if dark matter injects positrons into the Galactic medium, diffuse $\gamma$-rays offer the best channel for detection.

7 Conclusions

We have studied the prospects of indirectly detecting MeV DM by looking for sharp spectral features with a future MeV $\gamma$-ray mission.

We considered both prompt and secondary emission, for a large range of kinematically accessible two-body final states. We focused on emission from the inner Galaxy, where the effect of secondary $\gamma$-rays is especially important. Using DRAGON, we have consistently modeled for the first time all the relevant energy-loss and radiative processes relevant for leptonic final states at MeV energies. In addition, we considered the impact of diffusion and advection ($v=250\mathrm{\,km\,s^{-1}}$), with the latter dominating cosmic–ray transport. We found that secondary $\gamma$-ray emission, in particular in–flight annihilation and bremsstrahlung, will here often dominate the signal if $m_{\chi}\lesssim\mathcal{O}(100\mathrm{\,MeV})$ (see Fig. 5).

Projected upper limits for future MeV $\gamma$-ray missions, like the currently proposed e-ASTROGAM, were derived using a new approach based on Fisher forecasting with Poisson likelihoods. To this end, we assumed that background systematics for the future instruments will be at a similar level as for Fermi-LAT today.

For annihilation into leptonic final states, we found (Fig. 4) that diffuse $\gamma$-rays from the inner Galaxy can be used to probe the annihilation cross-section of sub-GeV DM at a level that is just a factor of a few weaker than current CMB constraints (for $s$-wave annihilation and standard halo profiles; for $p$-wave annihilation CMB constraints vanish and $\gamma$-rays are the stronger probe). Prompt $\gamma$-ray signals, such as monochromatic photons or box-like spectra, can be stringently tested throughout the entire MeV–GeV energy range, superseding current CMB limits by at least an order of magnitude (see Fig. 3). Remarkably, values of $\sigma v$ corresponding to thermal $p$–wave production in the early Universe can be probed below DM masses of $m_{\chi}\lesssim 20\mathrm{\,MeV}$.

In summary, future MeV missions such as proposed e-ASTROGAM have in general the potential to probe large parts of previously unexplored parameter space for annihilating sub-GeV DM, either improving existing constraints by at least one and up to several orders of magnitude, or detecting a sharp spectral feature possibly originating from DM annihilation.

Acknowledgments

We thank Gianfranco Bertone, Torsten Bringmann, Francesca Calore, Thomas Edwards, Mario Nicola Mazziotta, Christopher McCabe and Nicholas Rodd for useful discussions. SURFSara is thanked for use of the Lisa Compute Cluster. This research is funded by NWO through a GRAPPA-PhD fellowship (RB) and through an NWO VIDI research grant (CW).

Appendix A Annihilation spectra

Below we present analytic expression for the annihilation spectra of in-flight annihilation and final-state radiation.

A.1 In-flight annihilation

The differential annihilation cross section for a positron on an electron at rest is [see 79, 44]

[TABLE]

for $\gamma(1-\beta)\leq 2k-1\leq\gamma(1+\beta)$. Here $r_{e}$ is the classical radius of the electron, $m_{e}$ is the electron mass $\gamma$ is the Lorentz factor, $\beta=v/c$ and $k=E_{\gamma}/m_{e}$. From the kinematic bounds we can see that For non-relativistic positrons, i.e. $\gamma\rightarrow 1$ and $\beta\rightarrow 0$, photons from IfA are mono-energetic at $E_{\gamma}=m_{e}$. In the ultra-relativistic limit ($\gamma\rightarrow\infty$ and $\beta\rightarrow 1$) we get that $\gamma\left(1-\beta\right.)\rightarrow 0$. As such there is a lower limit to the photon energy of $m_{e}/2\leq E_{\gamma}$.

The production rate of photons from in-flight annihilation by a single positron is then,

[TABLE]

Here $N_{e^{-}}$ is the target density of electrons. Note that the differential cross section (Eq. A.1) is already weighted with a photon multiplicity of 2.

A.2 Final-state radiation

We present the expressions for final-state radiation applied in this paper.

Direct annihilation into $e^{+}e^{-}$.

In case of direct annihilation into fermions, the number of FSR photons per annihilation is: [62, 96, 64]:

[TABLE]

with $\sigma_{\mathrm{tot}}$ the total annihilation cross section into $e^{+}e^{-}$, $\alpha=1/137$ the fine-structure constant, $s=4m_{\chi}$ (roughly the center-of-mass energy for non-relativistic DM) and $s^{\prime}=4m_{\chi}(m_{\chi}-E_{\gamma})$.

Cascade annihilation

For a one-step cascade annihilation, $\chi\chi\rightarrow\phi\phi$, where subsequently $\phi\rightarrow e^{+}e^{-}$, the FSR spectrum is a little more complicated since the mediators, $\phi$, are boosted. An analytic expression can be derived by starting from Eq. A.3, but now for $\phi\rightarrow e^{+}e^{-}$, so $s=m_{\phi}^{2}$ and then boosting back to the rest frame of the DM. The resulting expression is [64],

[TABLE]

where $x=E_{\gamma}/m_{\chi}$ and $\mathrm{Li_{2}}(z)=\int_{z}^{0}\frac{\ln(1-t)}{t}dt$ is the dilogarithm.

Appendix B Cosmic-ray transport

In this appendix we discuss the results for different setups for the cosmic-ray transport. Corresponding limits for $\chi\chi\rightarrow e^{+}e^{-}$ are shown in Fig. 7. In addition, we also show how the secondary spectra change when the transport parameters are varied in Fig. 6. Our benchmark scenario includes an advective wind perpendicular to the plane with a velocity of $v_{\mathrm{wind}}=250\mathrm{\,km\,s^{-1}}$ and confined within the inner 3 kpc radius. The diffusion coefficient was modeled as a single power law. However, the behavior of the diffusion coefficient at low energies is unconstrained. We also study the case where there the diffusion coefficient flattens below 1.5 GeV, leading to faster diffusion at low energies. In this scenario, the diffusion coefficient as a function of momentum is modeled as a broken power law,

[TABLE]

where everything is as in Eq. 3.1, with the addition of the break momentum, $p_{b}=1.5\mathrm{\,GeV}$, and the the diffusion coefficient becomes energy independent below the break, i.e. $\delta_{l}=0$ (see Fig. 1). In addition, we consider scenarios in which there is no wind.

Figure 6 shows the spectra for a $m_{\chi}=60\mathrm{\,MeV}$ DM particle annihilating to an $e^{+}/e^{-}$ pair. Figure 6(a) is identical to Fig. 2 and depicts our benchmark model negligible diffusion at low momenta and an advective wind. In Fig. 6(b) the wind has been turned off. Naturally, final-state-radiation, being prompt, emission is unaffected. The secondary emission increases by a factor of a few compared to Fig. 6(a). Without advection and with single-power law diffusion the electrons and positrons loose their energy almost in-situ. Whereas in the presence of a wind they are likely to be advected out of our ROI before radiating. In the bottom two panels we show the results for the scenario in which the diffusion coefficient flattens below momenta of 1.5 GeV, again with (Fig. 6(c)) and without (Fig. 6(d)) wind. As can be seen, this has little impact on the result since in either the wind dominates transport as is the case in the bottom left plot, or cooling timescales are still faster than diffusion timescales (both bottom panels). All this is evident from Fig. 1

In addition, we compare projected upper limits the scenarios with and without wind and for single power-law diffusion in Fig. 7. We do not show limits for the scenario with a flat diffusion coefficient below 1.5 GeV as this only marginally affects the results, as discussed above. All assumptions on the instrument and background are identical to those in the main text. As mentioned above, in absence of the wind the secondary particles cool and radiate almost in-situ. Resulting in secondary emission being yet more important, with IfA being the dominant contribution to the limits upto $m_{\chi}\sim 100$ MeV.

Appendix C Spectra for different channels

Below we present example spectra for the different annihilation channels presented in this paper. All spectra are for benchmark setup and correspond to a window $|l|,|b|<5^{\circ}$. The annihilation cross section is everywhere set to $\left<\sigma v\right>=10^{-28}\mathrm{\,cm^{3}\,s^{-1}}$.

In Fig. 8(a) we show the spectrum from $\chi\chi\rightarrow e^{+}e^{-}$ for $m_{\chi}=60\mathrm{\,MeV}$. For the same dark matter mass we show the spectrum for the cascade channel in Fig. 8(b). The mediator mass is $m_{\phi}=5\mathrm{\,MeV}$. The FSR spectrum is now suppressed compared to the direct annihilation channel due to the smaller mass splitting between the mediator and and the electron, relative to the mass splitting between the electron and the dark matter particle. Also, the secondary spectra (IfA and bremsstrahlung) become more smooth since the injection spectrum of electrons and positrons gets spread out as a result of the boost of the mediator.

In the bottom panels we show the result for the muonic channel (Fig. 8(c)) at $m_{\chi}=110\mathrm{\,MeV}$. For the muonic channel FSR is the dominant component, which in this case is generated at two stages, when dark matter annihilates to a muon pair, and when the muon subsequently decays. Also, secondaries are more suppressed here compared to the direct or cascade annihilation channel, since the spectrum gets broadened as a result of the non-zero momentum of the muon and only one $e^{+}e^{-}$ pair is produced, compared to two in the cascade channel.

Finally, we show the spectrum for $\chi\chi\rightarrow\pi_{0}\gamma$ at $m_{\chi}=300\mathrm{\,MeV}$ in Fig. 8(d). The spectrum consists of a monochromatic line (which we gave a width of 2% for visual simplicity) and a box, since the pion is a scalar. The cross section is the same in all figures, but note the different vertical scale in Fig. 8(d) compared to the other figures.