section0em2.5em subsection2.5em2.5em

FERMILAB-PUB-21-457-T

# Reviving MeV-GeV Indirect Detection with Inelastic Dark Matter

Asher Berlina,b<EMAIL_ADDRESS>Gordan Krnjaica,c,d<EMAIL_ADDRESS>Elena
Pinettia,d<EMAIL_ADDRESS>aTheory Division, Fermi National Accelerator
Laboratory  bSuperconducting Quantum Materials and Systems Center (SQMS),
Fermi National Accelerator Laboratory  cUniversity of Chicago, Department of
Astronomy and Astrophysics  dUniversity of Chicago, Kavli Institute for
Cosmological Physics

###### Abstract

Thermal relic dark matter below $\sim 10\ \text{GeV}$ is excluded by cosmic
microwave background data if its annihilation to visible particles is
unsuppressed near the epoch of recombination. Usual model-building measures to
avoid this bound involve kinematically suppressing the annihilation rate in
the low-velocity limit, thereby yielding dim prospects for indirect detection
signatures at late times. In this work, we investigate a class of
cosmologically-viable sub-GeV thermal relics with late-time annihilation rates
that are detectable with existing and proposed telescopes across a wide range
of parameter space. We study a representative model of inelastic dark matter
featuring a stable state $\chi_{1}$ and a slightly heavier excited state
$\chi_{2}$ whose abundance is thermally depleted before recombination. Since
the kinetic energy of dark matter in the Milky Way is much larger than it is
during recombination, $\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$ upscattering can
efficiently regenerate a cosmologically long-lived Galactic population of
$\chi_{2}$, whose subsequent coannihilations with $\chi_{1}$ give rise to
observable gamma-rays in the $\sim 1\ \text{MeV}-100\ \text{MeV}$ energy
range. We find that proposed MeV gamma-ray telescopes, such as e-ASTROGAM,
AMEGO, and MAST, would be sensitive to much of the thermal relic parameter
space in this class of models and thereby enable both discovery and model
discrimination in the event of a signal at accelerator or direct detection
experiments.

##  I Introduction

Light dark matter (DM) has received much attention over the past decade as new
terrestrial experiments have been proposed to vastly expand laboratory
sensitivity to new physics below the electroweak scale – for reviews, see
Refs. Krnjaic  _et al._ (2022); Essig  _et al._ (2022). Amongst the various
possibilities, $\text{MeV}-\text{GeV}$ scale DM in the form of a light thermal
relic is especially motivated by basic principles of early-Universe
thermodynamics and known examples within the Standard Model (SM) Berlin  _et
al._ (2019).

Thermal DM is in equilibrium with the visible sector at early times and later
freezes out after becoming non-relativistic, such that number-changing
reactions occur less than once per Hubble time. Although such reactions remain
out of equilibrium, rare annihilation events can continue to inject energy
into the SM plasma. For instance, DM annihilations into electromagnetically-
charged particles can persist between recombination and reionization, yielding
observable distortions in the temperature anisotropies of the cosmic microwave
background (CMB). The lack of such deviations sharply constrains some of the
most predictive and well-motivated freeze-out models Padmanabhan and
Finkbeiner (2005); Slatyer (2016). This limit can be phrased in terms of the
thermally-averaged annihilation cross section $\langle\sigma v\rangle$ to
visible states, which at the time of recombination is constrained by the
Planck satellite at the level of Aghanim  _et al._ (2020)

$\langle\sigma v\rangle_{\text{cmb}}\lesssim 3\times 10^{-26}\
\text{cm}^{3}\,\text{s}^{-1}\ \bigg{(}\frac{m_{{}_{\text{DM}}}}{10\
\text{GeV}}\bigg{)}~{},$ (1)

where $m_{{}_{\text{DM}}}$ is the DM mass. Since standard thermal DM predicts
$\langle\sigma v\rangle\sim 10^{-26}\ \text{cm}^{3}\,\text{s}^{-1}$ at the
time of freeze-out when the DM was quasi-relativistic, models involving
annihilations to visible final states with velocity-independent ($s$-wave)
cross sections are excluded for masses below $\sim 10\ \text{GeV}$.

Of course, this cross section needs not be velocity-independent, and the CMB
bound is naturally evaded in models where the annihilation rate is suppressed
by the small DM velocity $v$ during recombination. For example, in models with
$p$-wave annihilation, $\langle\sigma v\rangle\propto v^{2}$, so that this
rate is significantly smaller at later times when DM is highly non-
relativistic. However, this approach111Although larger Galactic annihilation
rates can arise in other models of thermal relics involving resonant Feng and
Smolinsky (2017); Bernreuther  _et al._ (2021); Brahma  _et al._ (2023) or
forbidden annihilations Griest and Seckel (1991); D’Agnolo and Ruderman
(2015); D’Agnolo  _et al._ (2021), this is typically promising only when the
relative mass-splitting between various dark sector states is fixed to a
special value, corresponding to a strong tuning between mass scales that are
related to independent model parameters. also severely diminishes the
prospects for indirect detection in the Galaxy, since the local annihilation
rate is suppressed by the small virial velocity $v\sim 10^{-3}$, resulting in
$\sigma v<10^{-30}\ \text{cm}^{3}\,\text{s}^{-1}$ for standard thermal relics.
Since such small cross sections are beyond the reach of any near-future probe
(see Refs. Bartels  _et al._ (2017); Coogan  _et al._ (2021) and references
therein), there is strong motivation to identify well-motivated and testable
targets for indirect detection in the $\text{MeV}-\text{GeV}$ mass range.

Figure 1: A timeline of key events in the cosmic history of iDM. At early
times, coannihilations $\chi_{1}\chi_{2}\to\text{SM}\ \text{SM}$ freeze out at
a temperature much greater than the mass-splitting $T_{\text{fo}}\gg\delta$,
leaving behind roughly equal $\chi_{1,2}$ populations with a $\chi_{2}$ number
fraction $f_{2}=n_{2}/(n_{1}+n_{2})\simeq 1/2$. After freeze-out, $\chi_{2}$
continues to be depleted through DM-SM and DM-DM scattering. DM-SM scattering
additionally keeps the dark sector in kinetic equilibrium with the SM until
decoupling at a temperature of $T_{\text{kd}}$. At later times, the heavier
state continues to be depleted more rapidly through
$\chi_{2}\chi_{2}\to\chi_{1}\chi_{1}$ downscattering until decoupling at a DM
temperature of $T_{\chi}^{\text{dec}}$, such that the primordial fraction of
excited states is exponentially suppressed at the time of recombination. As a
result, CMB limits on DM annihilations are easily evaded. At much later times,
the virialization of DM in the Milky Way increases its velocity, such that
$\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$ upscattering partially restores the
abundance of long-lived heavier states and enhances the coannihilation rate
compared to its average cosmological value (red line). This increase implies a
local enhancement of $\chi_{1}\chi_{2}\to\text{SM}\ \text{SM}$ in the Milky
Way whenever the effective Galactic temperature of DM satisfies
$T_{\text{mw}}\gtrsim T_{\chi}^{\text{dec}}$. If the $\chi_{2}$ population
becomes appreciable in the Galaxy, then $\chi_{1}\chi_{2}\to\chi_{1}\chi_{2}$
elastic scattering can also play an important role in Galaxy evolution.

An alternative approach to evading the CMB limit on sub-GeV thermal DM
involves _changing_ the DM population between freeze-out and recombination, as
in models of “inelastic DM” (iDM) Tucker-Smith and Weiner (2001). In this
framework, such evolution occurs dynamically in a dark sector containing a
nearly-degenerate pseudo-Dirac DM pair, $\chi_{1}$ and $\chi_{2}$, with a
small mass-splitting and an off-diagonal coupling to the SM. The natural size
of the mass-splitting is near or below the overall DM mass scale if it arises
from the same dynamics responsible for generating the dark sector masses
Tucker-Smith and Weiner (2001); Arkani-Hamed  _et al._ (2009); Arkani-Hamed
and Weiner (2008). In this case, DM annihilation proceeds efficiently at early
times through coannihilations of $\chi_{1}$ and $\chi_{2}$ to SM final states,
but is exponentially suppressed before recombination as the Universe cools and
DM self-scattering continues to thermally deplete the abundance of the heavier
$\chi_{2}$ state in favor of the lighter $\chi_{1}$ state. A summary of the
cosmological history is depicted schematically in Fig. 1.

In this paper, we point out that the kinetic energy of $\chi_{1}$ in the Milky
Way is sufficiently large to undo this exponential suppression at late times
and thus regenerate a population of upscattered $\chi_{2}$ excited states that
are detectable through their subsequent coannihilation to visible states,
$\chi_{1}\chi_{2}\to\ \text{SM}\ \text{SM}$. In particular, we highlight the
existence of cosmologically-viable and theoretically well-motivated parameter
space involving mass-splittings of $\sim 1\ \text{eV}-100\ \text{eV}$, such
that the kinetic energy of $\text{MeV}-\text{GeV}$ scale DM in the Milky Way
is sufficient to partially or fully equilibrate a density of cosmologically
long-lived excited states. As a result, coannihilations that yield low-energy
gamma-rays occur at a rate comparable to (or a few orders of magnitude below)
the thermal value for $s$-wave annihilation, which is detectable with existing
and proposed indirect detection searches in the Milky Way. In contrast,
analogous signals of such processes are suppressed in the early Universe as
well as in less dense astrophysical bodies, such as dwarf satellite galaxies,
due to the smaller DM velocities involved.

We note that analogous ideas involving the regeneration of late-time indirect
detection signals have been previously explored within the context of prompt
decays of excited iDM states Finkbeiner and Weiner (2007); Chen  _et al._
(2010); Finkbeiner and Weiner (2016), late-time decays in a multi-component
hidden sector D’Eramo and Profumo (2018), and annihilations from asymmetric DM
oscillations Buckley and Profumo (2012). While such work is similar in spirit
to our study, these models typically involve additional particle content or
early-Universe dynamics compared to the model here. By contrast, we present a
simple scenario with standard cosmological assumptions where the same
interaction that sets the DM thermal relic abundance is also responsible for
indirect detection signals in the Galaxy.

The remainder of this paper is organized as follows. In Secs. II and III, we
provide an overview of iDM models and their cosmological history. In Sec.
III.4, these results are used to determine the primordial fraction of excited
states at late times. In Sec. IV, we study the Galactic dynamics of this
scenario and calculate the excited DM population in the Milky Way. These
results are then used in Secs. V and VI to derive limits from existing
indirect detection searches as well as prospects for observing such a
population in low-energy gamma-rays with proposed satellites. Finally, we
offer concluding remarks in Sec. VII and discuss directions for future work. A
series of appendices, referred to throughout this work, provides additional
technical details.

##  II Model Overview

In this section, we outline a concrete model of iDM. The DM content of this
model involves a Dirac pair of two-component Weyl fermions, $\psi_{1}$ and
$\psi_{2}$, oppositely charged under a new spontaneously broken
$U(1)^{\prime}$ symmetry with gauge interaction

$\mathscr{L}\supset
e^{\prime}A^{\prime}_{\mu}\big{(}\psi_{1}^{\dagger}\bar{\sigma}^{\mu}\psi_{1}-\psi_{2}^{\dagger}\bar{\sigma}^{\mu}\psi_{2}\big{)}~{},$
(2)

where the dark photon $A^{\prime}$ is the $U(1)^{\prime}$ gauge boson with
mass $m_{A^{\prime}}$ and gauge coupling $e^{\prime}$. The fermion masses of
this theory arise from the Lagrangian terms

$\mathscr{L}\supset-
m_{D}\,\psi_{1}\,\psi_{2}-\frac{1}{2}\big{(}\delta_{1}\,\psi_{1}\psi_{1}+\,\delta_{2}\,\psi_{2}\psi_{2}\big{)}+\text{h.c.}~{},$
(3)

where $m_{D}$ is a gauge-invariant Dirac mass and $\delta_{1,2}$ are
$U(1)^{\prime}$ breaking Majorana masses naturally222Indeed, the Majorana
terms are small if they arise from higher-dimensional couplings to the dark
Higgs $H^{\prime}$ responsible for spontaneous breaking in the dark sector
(e.g., couplings of the form $H^{\prime\,2}\,\psi_{1,2}^{2}/\Lambda$ where
$\Lambda$ is governed by the mass scale of a heavy singlet that mixes with
$H^{\prime}$) Tucker-Smith and Weiner (2001); Arkani-Hamed  _et al._ (2009);
Arkani-Hamed and Weiner (2008). More generally, $\delta_{1,2}\ll m_{D}$ is
technically natural as the theory enjoys an enhanced global symmetry
(analogous to SM lepton number) when Majorana masses are absent. satisfying
$\delta_{1,2}\ll m_{D}$. After diagonalizing this system, the mass eigenstates
correspond to a pseudo-Dirac pair,

$\chi_{1}\simeq\frac{i}{\sqrt{2}}\,(\psi_{1}-\psi_{2})~{}~{},~{}~{}\chi_{2}\simeq\frac{1}{\sqrt{2}}\,(\psi_{1}+\psi_{2})~{},$
(4)

with nearly degenerate masses

$m_{1,2}\simeq m_{D}\mp\frac{1}{2}(\delta_{1}+\delta_{2})$ (5)

split by the small amount $\delta\equiv\delta_{1}+\delta_{2}\ll m_{1,2}$. In
terms of mass eigenstates, the leading interaction of Eq. 2 is off-diagonal,

$\mathscr{L}\supset
ie^{\prime}A^{\prime}_{\mu}\,\bar{\chi}_{1}\gamma^{\mu}\chi_{2}+\mathcal{O}\bigg{(}\frac{\delta_{1,2}}{m_{D}}\bigg{)}~{},$
(6)

where we have written $\chi_{1,2}$ as four-component Majorana spinors.

Generically, $\chi_{1}$ and $\chi_{2}$ also couple diagonally to the
$A^{\prime}$, but these interactions are suppressed by $\delta/m_{1}\ll 1$ and
exactly vanish if $\delta_{1}=\delta_{2}$ due to the enhanced charge
conjugation symmetry $\chi_{1,2}\to\mp\chi_{1,2}$, $A^{\prime}\to-A^{\prime}$.
Without loss of essential generality, throughout this paper we assume that
$\chi_{1,2}$ couple purely off-diagonally to the dark photon, as in Eq. 6,
which governs the dynamics within the dark sector. For instance, Eq. 6
facilitates upscattering $\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$ as well as
_elastic_ self-scattering $\chi_{1}\chi_{2}\to\chi_{1}\chi_{2}$ between DM
particles, both of which can have significant impacts on the Galactic
population, as we will discuss below.

The dark photon coupling to the electromagnetic current arises from a small
kinetic mixing between $A^{\prime}$ and the SM photon field $A$,

$\mathscr{L}\supset\frac{\epsilon}{2}\,F_{\mu\nu}^{\prime}\,F^{\mu\nu}+A_{\mu}J_{\text{em}}^{\mu}~{},$
(7)

where $F^{\prime}$ and $F$ are the dark and visible field strengths,
respectively, $J_{\text{em}}$ is the SM electromagnetic current, and
$\epsilon\ll 1$ is the dimensionless parameter governing the strength of
$A^{\prime}$-$A$ mixing. After rotating away this mixing via $A\to A+\epsilon
A^{\prime}$, SM fields acquire a small coupling to the dark photon through the
induced interaction $\epsilon\,A^{\prime}_{\mu}J^{\mu}_{\text{em}}$ Holdom
(1986).

##  III Cosmological History

In this section, we review the cosmological history of sub-GeV iDM. The
discussion is organized into different subsections, ordered chronologically
according to when various processes decouple in the early Universe in the
parameter space of interest. In particular, we discuss $\text{DM}-\text{SM}$
chemical decoupling, $\text{DM}-\text{SM}$ kinetic decoupling, and
$\chi_{1}-\chi_{2}$ chemical decoupling, respectively, as depicted in Fig. 1.
For previous work investigating such cosmological scenarios, see, e.g., Refs.
Baryakhtar  _et al._ (2022); Carrillo González and Toro (2022).

###  III.1 $\text{DM}-\text{SM}$ Chemical Decoupling

Here we discuss DM thermal freeze-out, the process of $\text{DM}-\text{SM}$
chemical decoupling which sets the total $\chi_{1,2}$ density in the early
Universe. If the kinetic mixing parameter satisfies

$\epsilon\gtrsim 10^{-8}\,\sqrt{\frac{m_{A^{\prime}}}{\text{GeV}}}~{}~{},$ (8)

then $A^{\prime}\leftrightarrow\text{SM}\ \text{SM}$ decays and inverse-decays
equilibrate the two sectors, seeding an initially large thermal density of
dark sector states Evans  _et al._ (2018).

A priori, there is no preferred ordering of DM and dark photon masses.
However, if $m_{A^{\prime}}<m_{1}$ and $e^{\prime}\gg\epsilon$, then DM
freeze-out proceeds predominantly through $\chi_{1}\chi_{1}\to
A^{\prime}A^{\prime}$ annihilation followed by $A^{\prime}$ decays to visible
SM particles. Since this process is $s$-wave, the CMB limit from Eq. 1
robustly excludes this mass ordering in the sub-GeV mass range.

By contrast, if $m_{A^{\prime}}>m_{1}+m_{2}$, then DM annihilations to on-
shell dark photons are kinematically forbidden, and instead the DM relic
density is governed by the coannihilation $\chi_{1}\chi_{2}\to\text{SM}\
\text{SM}$, as shown in Fig. 1. In the non-relativistic limit, the cross
section for this process is Berlin  _et al._ (2014, 2019)

$\langle\sigma
v\rangle_{\text{ann}}\simeq\frac{16\pi\alpha\alpha^{\prime}\epsilon^{2}m_{1}^{2}}{(4m_{1}^{2}-m_{A^{\prime}}^{2})^{2}}~{},$
(9)

where $\alpha$ and $\alpha^{\prime}=e^{\prime\,2}/4\pi$ are the
electromagnetic and $U(1)^{\prime}$ fine-structure constants, respectively. In
our calculations, we incorporate all kinematically accessible channels,
including hadronic final states, as described in Refs. Izaguirre  _et al._
(2016); Berlin and Kling (2019). Although this cross section is unsuppressed
in the low-velocity limit, the CMB bound is alleviated because the $\chi_{2}$
abundance is thermally depleted before recombination for sufficiently large
splittings, thereby shutting off $\chi_{1}\chi_{2}\to\ \text{SM}\ \text{SM}$
at late times. This is discussed in more detail below.

When the temperature of the early Universe cools to $T=T_{\text{fo}}\sim
m_{1}/10$, the total comoving DM density freezes out to a fixed value. Thus,
for $\delta\ll m_{1}/10$, the mass-splitting can be neglected throughout this
process, and we follow the calculations of Ref. Berlin  _et al._ (2019) to
determine the initial relic density of $\chi_{1,2}$ particles. In particular,
for $\delta\ll m_{1,2}\lesssim 100\ \text{MeV}$ and
$m_{A^{\prime}}\gtrsim\text{few}\times m_{1}$, the $\chi_{1,2}$ abundance at
freeze-out satisfies

$f_{\chi}\simeq\bigg{(}\frac{0.5}{\alpha^{\prime}}\bigg{)}\bigg{(}\frac{10^{-5}}{\epsilon}\bigg{)}^{2}\bigg{(}\frac{m_{1}}{100\
\text{MeV}}\bigg{)}^{2}\bigg{(}\frac{m_{A^{\prime}}}{3\,m_{1}}\bigg{)}^{4}~{},$
(10)

where $f_{\chi}$ is the fractional abundance of DM that $\chi_{1}$ and
$\chi_{2}$ jointly constitute.

Throughout this work, we focus on the parameter space that satisfies the
assumptions leading to Eq. 10, fixing to the representative values
$m_{A^{\prime}}/m_{1}=3$ and $\alpha^{\prime}=0.5$, and assuming a standard
freeze-out cosmology for $\chi_{1}$ and $\chi_{2}$.333Thermal freeze-out
models with $m_{A^{\prime}}/m_{1}\lesssim 1$ are typically excluded by
measurements of the CMB, as mentioned above, whereas scenarios with
$m_{A^{\prime}}/m_{1}\gg 3$ or $\alpha^{\prime}\ll 1$ are strongly constrained
by accelerator searches Berlin  _et al._ (2019, 2020). We also fix the value
of $\epsilon$ by imposing $f_{\chi}=1$ such that $\chi_{1,2}$ account for the
entire observed DM abundance, unless mentioned otherwise (in Sec. VI we will
also investigate DM subcomponents, corresponding to $f_{\chi}<1$). In later
sections, we vary the iDM mass-splitting $\delta$ to explore a wide range of
possible late-time implications of the $\chi_{2}$ population. However, in the
parameter space of interest, it is always the case that $\delta\ll m_{1}$, so
that Eq. 10 is unaffected.

Furthermore, since we always work in the regime where
$\delta<m_{A^{\prime}},2m_{e}$, the decays $\chi_{2}\to\chi_{1}A^{\prime}$ and
$\chi_{2}\to\chi_{1}e^{+}e^{-}$ are kinematically forbidden and $\chi_{2}$ is
stable on cosmological timescales; we neglect the highly subdominant
$\chi_{2}\to\chi_{1}+3\gamma$ and $\chi_{2}\to\chi_{1}\nu\bar{\nu}$ decay
channels, since the corresponding lifetimes are exponentially longer than the
age of the Universe Batell  _et al._ (2009).

###  III.2 $\text{DM}-\text{SM}$ Kinetic Decoupling

After $\text{DM}-\text{SM}$ chemical decoupling, $\text{DM}-\text{SM}$
scattering $\chi_{1}\,\ell^{\pm}\leftrightarrow\chi_{2}\,\ell^{\pm}$ remains
in equilibrium due to the large thermal density of SM leptons,444Note that
there are also similar interactions with quarks and hadrons, but since their
densities are exponentially suppressed at temperatures $T\lesssim 100\
\text{MeV}$, such interactions are subdominant.  as shown in Fig. 1. These
processes do not change the total DM number, but as long as they are in
equilibrium, the DM temperature, $T_{\chi}$, tracks that of the SM, $T$, such
that both approximately evolve inversely with the cosmic scale factor
$T_{\chi}=T\propto a^{-1}$.

Since our focus is on sub-GeV DM, the dominant process that maintains kinetic
equilibrium after freeze-out is inelastic scattering off of electrons and
positrons, whose abundance vastly exceeds that of any other charged species at
temperatures $T\lesssim 100\ \text{MeV}$. The full thermally-averaged
thermalization rate for $\chi_{1}\,e\leftrightarrow\chi_{2}\,e$ is derived in
Appendix A, but here we briefly summarize its limiting forms. At both high and
low temperatures compared to the electron mass, this rate is well approximated
by

$\Gamma_{\text{kd}}\simeq\frac{8\alpha\alpha^{\prime}\epsilon^{2}T^{2}}{3m_{A^{\prime}}^{4}m_{1}}\times\begin{cases}\frac{31}{63}\,\pi^{5}\,T^{4}&(T\gg
m_{e})\\\ \frac{16}{\pi}\,m_{e}^{4}\,e^{-m_{e}/T}&(T\ll m_{e})~{}.\end{cases}$
(11)

The DM and SM sectors kinetically decouple from each other at a temperature of
$T_{\chi}=T=T_{\text{kd}}$ when the typical DM particle no longer efficiently
transfers momentum with an electron in a Hubble time. This temperature can be
approximated by solving $\Gamma_{\text{kd}}(T_{\text{kd}})=H(T_{\text{kd}})$,
where $H$ is the Hubble expansion rate. After this point, the DM temperature
no longer tracks that of the SM plasma, instead redshifting more rapidly as an
isolated non-relativistic sector with a temperature scaling of
$T_{\chi}\propto a^{-2}$.

For $m_{1}\lesssim 1\ \text{GeV}$ and our choice of representative model
parameters in Eq. 10, $\Gamma_{\text{kd}}\gg H$ at temperatures comparable to
or greater than the electron mass. Thus,
$\chi_{1}e^{\pm}\leftrightarrow\chi_{2}e^{\pm}$ scattering steadily depletes
the $\chi_{2}$ population until $T\sim m_{e}$ when the electron/positron
number density becomes Boltzmann suppressed and the scattering rate falls out
of equilibrium. Note that for $m_{1}\lesssim 10\ \text{MeV}$, DM
annihilations, as discussed in Sec. III.1, decouple after
$\text{DM}-\text{SM}$ scattering, so that $T_{\text{kd}}>T_{\text{fo}}$ and
the two earliest processes in Fig. 1 decouple in the opposite order relative
to what is shown. In this case, once DM annihilations to the SM freeze out,
the temperatures of the two sectors evolve independently. Generally, we expect
this to only hold approximately, and we reserve a more careful numerical
treatment, involving solving the relevant Boltzmann equations, to future work.

###  III.3 $\text{DM}-\text{DM}$ Chemical Decoupling

DM self-scattering $\chi_{1}\chi_{1}\leftrightarrow\chi_{2}\chi_{2}$ (as shown
in Fig. 1) is unsuppressed by the small coupling $\epsilon$, and thus can
remain in equilibrium long after the processes of the previous two subsections
have already decoupled. This implies that even for splittings as small as
$\delta\sim\text{eV}$, $\chi_{2}$ can be thermally depleted well before
recombination once $T_{\chi}\lesssim\delta$.

As derived in Appendix B, the thermally-averaged rate for
$\chi_{2}\,\chi_{2}\to\chi_{1}\,\chi_{1}$ scattering is given by

$\Gamma_{\chi}\simeq\frac{8\pi\,\alpha^{\prime\,2}\,m_{1}^{3/2}}{m_{A^{\prime}}^{4}}\,n_{1}\,e^{-\delta/T_{\chi}}\,\max\bigg{(}\frac{\delta}{2}\,,\,\frac{T_{\chi}}{\pi}\bigg{)}^{1/2}~{},$
(12)

where $n_{1}$ is the $\chi_{1}$ number density.555We note that Eq. 12 agrees
with the expressions in Ref. Carrillo González and Toro (2022), aside from a
factor of $\pi$ discrepancy for $T_{\chi}\lesssim\delta$. We estimate the DM
temperature $T_{\chi}^{\text{dec}}$ at which $\chi_{1}$ and $\chi_{2}$
chemically decouple from each other by evaluating
$\Gamma_{\chi}(T_{\chi}^{\text{dec}})=H$. After decoupling, the primordial
abundance of excited states is fixed to

$f_{2}=\frac{n_{2}}{n_{1}+n_{2}}\simeq\frac{e^{-\delta/T_{\chi}^{\text{dec}}}}{1+e^{-\delta/T_{\chi}^{\text{dec}}}}~{},$
(13)

where we have used detailed balance to relate the number densities of
$\chi_{2}$ and $\chi_{1}$ just prior to decoupling, $n_{2}\simeq
e^{-\delta/T_{\chi}}\,n_{1}$.

###  III.4 Primordial Excited Fraction

Figure 2: Contours of fractional $\chi_{2}$ abundance
$f_{2}=n_{2}/(n_{1}+n_{2})$ set by $\chi_{1}-\chi_{2}$ decoupling in the early
Universe, as a function of $m_{1}$ and $\delta$. Here, we have fixed
$m_{A^{\prime}}/m_{1}=3$, $\alpha^{\prime}=0.5$, and $\epsilon$ such that the
freeze-out abundance of $\chi_{1,2}$ agrees with the observed DM density. The
gray shaded region corresponds to parameter space excluded by CMB annihilation
bounds Aghanim  _et al._ (2020).

Using the formalism of Secs. III.1, III.2, and III.3, we determine the total
$\chi_{1,2}$ density, temperature evolution, and $\chi_{2}$ fractional density
$f_{2}$, respectively. In Fig. 2, we show $f_{2}$ as a function of $m_{1}$ and
$\delta$. For fixed $\delta$, heavier dark sector masses suppress the rate
$\Gamma_{\chi}$ of Eq. 12, such that $\chi_{1}-\chi_{2}$ interconversion
decouples at earlier times, enhancing the relative $\chi_{2}$ density. On the
other hand, for larger splittings, the relative abundance of excited states is
increasingly Boltzmann suppressed, yielding smaller values of $f_{2}$. This
cosmological value of $f_{2}$ is fixed at late times, as shown schematically
in Fig. 1 by the lower plateau after $T_{\chi}=T_{\chi}^{\text{dec}}$.

This primordial $\chi_{2}$ fraction governs the rate for
$\chi_{1}\chi_{2}\to\text{SM}\ \text{SM}$ near the time of recombination,
which is strongly constrained by CMB data. In particular, Planck bounds such
coannihilations to $e^{+}e^{-}$ final states according to Aghanim  _et al._
(2020)

$2f_{\chi}^{2}f_{1}f_{2}\langle\sigma v\rangle_{\text{ann}}\lesssim 2\times
10^{-26}\
\text{cm}^{3}\,\text{s}^{-1}\,\bigg{(}\frac{m_{1}}{30\,\text{GeV}}\bigg{)}~{},$
(14)

where the factor of two on the left-hand-side of Eq. 14 accounts for the fact
that the limit in Ref. Aghanim  _et al._ (2020) assumed annihilation of self-
conjugate DM particles, whereas coannihilations of $\chi_{1}$ and $\chi_{2}$
involve distinct species. We have also rescaled the limit by
$f_{\chi}^{2}f_{1}f_{2}$, where $f_{\chi}$ from Eq. 10 accounts for the DM
fraction that $\chi_{1}$ and $\chi_{2}$ jointly constitute, and
$f_{1}=1-f_{2}$ is the relative $\chi_{1}$ fraction defined in analogy with
$f_{2}$ in Eq. 13. Annihilations to muon, tau, and hadronic final states are
also similarly constrained, but are weaker by roughly a factor of $\sim 1/3$
due to their smaller efficiency in transferring energy into the intergalactic
medium near the time of recombination Slatyer (2016). The resulting constraint
from the CMB is shown in Fig. 2 as the shaded gray region.

##  IV Galactic Population

After dark sector scattering decouples at $T_{\chi}^{\text{dec}}$, the average
cosmological fraction $f_{2}$ of the heavier $\chi_{2}$ state remains fixed
for all subsequent times as the DM temperature continues to fall as $a^{-2}$
due to Hubble expansion. However, in overdense regions where galactic
structure forms, the local DM temperature can vastly exceed its average
cosmological value due to virialization. Thus, in galaxies and clusters where
the local temperature satisfies $T_{\chi}\gtrsim\delta$, there is an
additional late-time population of excited states generated by
$\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$ upscattering, which can exceed the
cosmological fraction $f_{2}$ computed in the previous section. In this
section, we discuss how such a population can arise in the Milky Way.

###  IV.1 Galactic Model

Since the Galactic $\chi_{2}$ fraction is exponentially sensitive to the DM
velocity, we begin by describing a toy model for the Milky Way. We assume that
upon galaxy formation, the $\chi_{2}$ population has the same cosmological
relic fraction $f_{2}$ set by dynamics in the early Universe (see Sec. III)
and starts off with the same type of Galactic density profile as the
$\chi_{1}$ population.

To simplify our treatment, we approximate the phase-space of DM in the Milky
Way as following that of a Maxwell-Boltzmann distribution, with an _effective_
temperature $T_{\text{mw}}$ that is a function of the galactocentric radius
$r$. In particular, using the virial theorem, we take

$T_{\text{mw}}(r)=\frac{GM_{\text{enc}}(r)}{3r}\,m_{1}~{},$ (15)

where $M_{\text{enc}}(r)$ is the total mass enclosed within $r$.
$M_{\text{enc}}$ is determined directly from our assumed DM density profile as
well as the baryonic mass density $\rho_{B}$. We approximate the stellar
bulge, thin disk, and thick disk contributions to $\rho_{B}$ from the best-fit
spherically-symmetric Hernquist profile model advocated for in Refs. McMillan
(2011); Kaplinghat  _et al._ (2014),

$\rho_{B}(r)=\frac{\rho_{B0}\,r_{0}^{4}}{r(r+r_{0})^{3}}~{},$ (16)

where $\rho_{B0}=26\ \text{GeV}\,\text{cm}^{-3}$ and $r_{0}=2.7\ \text{kpc}$.

Later in Sec. V, we will investigate the impact of various choices for the DM
density profile. However, for concreteness, in this section we model the
initial $\chi_{1}$ mass density with a Navarro-Frenk-White (NFW) profile
Navarro  _et al._ (1997)

$\rho_{1}(r)=\frac{\rho_{s}}{(r/r_{s})\,(1+r/r_{s})^{2}}~{},$ (17)

where $r_{s}=20\ \text{kpc}$ is the scale radius and $\rho_{s}$ is normalized
to yield a local DM density of $0.4\ \text{GeV}\,\text{cm}^{-3}$ at the
position of the Solar System, $r=r_{\odot}\simeq 8.5\ \text{kpc}$.

###  IV.2 Galactic Evolution

The time-evolution of the Galactic $\chi_{2}$ population is governed by the
Boltzmann equation

$\frac{\partial{\cal F}_{2}}{\partial t}+\vec{v}\cdot\vec{\nabla}{\cal
F}_{2}+\vec{g}\cdot\frac{\partial{\cal F}_{2}}{\partial\vec{v}}=C[{\cal
F}_{1,2}]~{},$ (18)

where ${\cal F}_{1,2}$ are the phase space densities for $\chi_{1,2}$,
$\vec{v}$ is the $\chi_{2}$ velocity, $\vec{g}$ is the Galactic gravitational
acceleration field, and the collision term $C[{\cal F}_{1,2}]$ accounts for
all $\chi_{1,2}$ scattering processes. Since solving Eq. 18 requires detailed
numerical simulations, we simplify the analysis by factorizing the processes
of _inelastic_ $\chi_{2}$ production and gravitational evolution, examining
upscattering with a simple toy model in this section. Modifications to the
$\chi_{1,2}$ distributions due to gravitational dynamics and elastic
scattering will be discussed later in Sec. V.2, where we consider a wide range
of possible scenarios.

Figure 3: Radial regions in the Milky Way in which inelastic dark scattering
$\chi_{1}\chi_{1}\leftrightarrow\chi_{2}\chi_{2}$ reenters equilibrium over
the age of the Galaxy $t_{\text{mw}}\simeq 13.5\ \text{Gyr}$, as a function of
$m_{1}$ and for various representative choices of $\delta$, fixing
$m_{A^{\prime}}/m_{1}=3$ and $\alpha^{\prime}=0.5$. We have also assumed an
initial NFW profile for $\chi_{1}$. For each choice of $m_{1}$ and $\delta$,
there is a shell of galactocentric radii within which equilibration occurs.

We approximate $\chi_{2}$ production by taking the zeroth moment of Eq. 18 and
neglecting gradient and gravitational terms, which amounts to solving

$\frac{\partial n_{2}}{\partial t}=n_{1}^{2}\,\langle\sigma v\rangle_{1\to
2}-n_{2}^{2}\,\langle\sigma v\rangle_{2\to 1}~{},$ (19)

where $n_{1,2}$ is the Galactic number density of $\chi_{1,2}$. Here, the
thermally-averaged cross section for $\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$
upscattering is

$\langle\sigma v\rangle_{1\to 2}=e^{-2\delta/T_{\text{mw}}}\,\langle\sigma
v\rangle_{2\to 1}~{},$ (20)

where $\langle\sigma v\rangle_{2\to 1}$ is the thermally-averaged cross
section for $\chi_{2}\chi_{2}\to\chi_{1}\chi_{1}$ downscattering. As derived
in Appendix B, $\langle\sigma v\rangle_{2\to 1}$ is approximately

$\langle\sigma v\rangle_{2\to
1}\simeq\frac{8\pi\,\alpha^{\prime\,2}\,m_{1}^{3/2}}{m_{A^{\prime}}^{4}}\,\,\max\bigg{(}\frac{\delta}{2}\,,\,\frac{T_{\text{mw}}}{\pi}\bigg{)}^{1/2}~{},$
(21)

From Eq. 19, $n_{2}$ asymptotes to a constant value once it evolves to
$n_{2}\simeq e^{-\delta/T_{\text{mw}}}\,n_{1}$, as required by detailed
balance. In the limit that the initial $\chi_{2}$ density is small, we can
therefore approximate the solution of Eq. 19 as

$n_{2}\simeq
n_{1}\,\min\big{(}N_{\text{scatt}}\,,\,e^{-\delta/T_{\text{mw}}}\big{)}~{},$
(22)

where $N_{\text{scatt}}=n_{1}\,\langle\sigma v\rangle_{1\to 2}\,t$ is the
number of upscatters per $\chi_{1}$ particle after time $t$. Note that
$\chi_{2}$ approaches detailed balance once the average number of upscatters
is $N_{\text{scatt}}\sim e^{-\delta/T_{\text{mw}}}$, which occurs on a
characteristic timescale of

$t\sim\frac{e^{\delta/T_{\text{mw}}}}{n_{1}\,\langle\sigma v\rangle_{2\to
1}}~{}.$ (23)

Due to its strong temperature dependence, the upscattering cross section
$\langle\sigma v\rangle_{1\to 2}$ is peaked around a characteristic radial
range in the Galaxy, driven by the fact that $T_{\text{mw}}$ peaks at
intermediate radii. In Fig. 3, we show the radial regions of the Milky Way in
which $\chi_{2}$ reenters chemical equilibrium with $\chi_{1}$ over the age of
the Galaxy, $t_{\text{mw}}\simeq 13.5\ \text{Gyr}$, as a function of $m_{1}$
and for various choices of the splitting $\delta$. We see that for mass-
splittings $\delta\lesssim 10\ \text{eV}$, $\chi_{2}$ is able to reenter
chemical equilibrium over a significant portion of the Galaxy.

We now calculate the Galactic density profile of $\chi_{1}$ and $\chi_{2}$ at
$t=t_{\text{mw}}$ by numerically solving Eq. 19, along with
$n_{1}+n_{2}=\text{constant}$, at each radius $r$, starting from the initial
condition at $t=0$ where

$n_{1,2}(0)=f_{1,2}\,f_{\chi}\,\rho_{{}_{\text{DM}}}/m_{1,2}~{},$ (24)

for a given Galactic DM profile $\rho_{{}_{\text{DM}}}$. Note that this
treatment neglects the impact of additional effects, including elastic
$\chi_{1}\chi_{2}\to\chi_{1}\chi_{2}$ scattering and gravitational dynamics,
which will be addressed later in Sec. V.2. In Fig. 4, we show the resulting
Galactic mass density profile of $\chi_{2}$, $\rho_{2}=m_{2}\,n_{2}$, fixing
$m_{1}=50\ \text{MeV}$ and $\delta=5\ \text{eV}$ and for two different
assumptions regarding the Galactic evolution of excited states. In particular,
for the solid line in Fig. 4, we take $\chi_{2}$ to be produced along orbits
of fixed radius. By contrast, for the dashed line we instead assume that after
production, the entire Galactic $\chi_{2}$ population readjusts to an NFW
profile with normalization fixed to conserve the total $\chi_{2}$ number.
These two scenarios are referred to as “unvirialized” and “virialized,”
respectively.666Note that our use of “virialization” is simply meant to imply
the process of gravitational relaxation/evolution. For reference, we also show
as a dotted line the initial primordial $\chi_{2}$ profile (see Eq. 24) and as
a dot-dashed line a standard NFW DM density profile (scaled down by a factor
of two). The peak-like feature near $r\sim 1\ \text{kpc}$ shows the growth to
the $\chi_{2}$ population generated from $\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$
Galactic upscattering at late times. Note that for much larger radii, the
distribution does not significantly evolve from its original primordial value
since the low density of $\chi_{1}$ particles suppresses the likelihood to
upscatter over the age of the Galaxy. Instead, for much smaller radii, the
small Galactic temperature favors downscattering compared to upscattering,
reducing the density of excited states in that region.

Figure 4: Radial mass density profiles of $\chi_{2}$ in the Milky Way, for
$m_{1}=50\ \text{MeV}$, $\delta=5\ \text{eV}$, $m_{A^{\prime}}/m_{1}=3$,
$\alpha^{\prime}=0.5$, and assuming an initial NFW profile for $\chi_{1}$ and
$\chi_{2}$. For the solid line labelled “unvirialized,” we assume that
$\chi_{2}$ remains at the same radius as the progenitor $\chi_{1}$ particle.
For the dashed line labelled “revirialized,” we instead assume that the
$\chi_{2}$ population readjusts to an NFW profile with normalization fixed to
conserve particle number. For reference, we also show as a dotted line the
initial primordial $\chi_{2}$ profile (see Eq. 24) and as a dot-dashed line a
standard NFW DM density profile, scaled down by a factor of two.

##  V Indirect Detection

The last section discussed the Galactic density profiles for $\chi_{1}$ and
$\chi_{2}$. Here, we use these results to determine the resulting flux of
visible particles from DM coannihilations $\chi_{1}\chi_{2}\to\text{SM}\
\text{SM}$. We begin by providing some general formalism in Sec. V.1. In Sec.
V.2, we then investigate the sensitivity of our final results to the assumed
form of the $\chi_{1}$ and $\chi_{2}$ density profiles.

###  V.1 General Formalism

In conventional indirect detection studies, the flux $\Phi_{\gamma}$ of
photons from DM annihilations is determined by

$\frac{d\Phi_{\gamma}}{dE_{\gamma}}=c_{\chi}\,\frac{\langle\sigma
v\rangle_{\text{ann}}}{4\pi\,m_{{}_{\text{DM}}}^{2}}\,\frac{dN_{\gamma}}{dE_{\gamma}}\,J~{},$
(25)

where $\langle\sigma v\rangle_{\text{ann}}$ is the thermally-averaged
annihilation cross section in the Galaxy, $dN_{\gamma}/dE_{\gamma}$ is the
number of photons per annihilation event per unit energy $E_{\gamma}$, and the
$J$-factor depends on the Galactic mass density of DM,
$\rho_{{}_{\text{DM}}}$. The constant $c_{\chi}$ incorporates additional
factors arising from the counting of DM species. For example, $c_{\chi}=1/2$
for Majorana DM in order to not overcount distinct pairs of self-conjugate
particles, and $c_{\chi}=1/4$ for Dirac DM, since annihilations involve
interactions between distinct particles and antiparticles, each of which
account for half of the total DM density. For standard DM annihilations, the
$J$-factor is given by

$J=\int d\Omega\,d\ell~{}\rho_{{}_{\text{DM}}}(r)^{2}~{},$ (26)

where $\Omega$ and $\ell$ are solid angle and radial coordinates with respect
to a local observer, which are related to the galactocentric radial coordinate
by

$r^{2}=r_{\odot}^{2}+\ell^{2}-2\,r_{\odot}\,\ell\cos{\theta}~{}.$ (27)

In our iDM scenario, there are two distinct DM species that coannihilate to
yield a flux of SM particles, similar to Dirac DM. However, unlike the Dirac
case, the $\chi_{2}$ population that arises from Galactic upscattering can
have a distinct profile, such that we define the $J$-factor for iDM as

$J_{\text{iDM}}\equiv\int d\Omega\,d\ell~{}\rho_{1}(r)\,\rho_{2}(r)~{},$ (28)

where the mass density profiles for $\chi_{1}$ and $\chi_{2}$ are determined
from Sec. IV.2. Note that in the limit where $m_{1}\simeq m_{2}$ and
$\rho_{1}(r)\simeq\rho_{2}(r)\simeq\rho_{{}_{\text{DM}}}/2$,
$J_{\text{iDM}}\simeq J/4$. With this convention for $J_{\text{iDM}}$,
$c_{\chi}=1$ for iDM in Eq. 25.

Most indirect detection studies of the Milky Way estimate the reach of various
existing or future telescopes by presenting the sensitivity to the
annihilation cross section $\langle\sigma v\rangle$, assuming an NFW profile
of self-conjugate DM particles. To make direct contact with these studies, we
define an _effective_ iDM annihilation cross section,

$\langle\sigma v\rangle_{\text{eff}}\equiv 2\,\langle\sigma
v\rangle_{\text{ann}}~{}\bigg{(}\frac{J_{\text{iDM}}}{J_{\text{NFW}}}\bigg{)}~{},$
(29)

where $J_{\text{NFW}}$ is defined as in Eq. 26 for an NFW DM profile. The
numerical prefactor in Eq. 29 is fixed so that a telescope that is sensitive
to self-conjugate (e.g., Majorana) DM annihilations with an NFW density
profile and an annihilation cross section of $\langle\sigma
v\rangle_{\text{eff}}$ is also sensitive to iDM particles of the same mass and
annihilation channels.

###  V.2 $J$-factor Profile-Dependence

Figure 5: As solid lines, the Milky Way $J$-factor for iDM (see Eq. 28)
normalized by that of conventional NFW-distributed DM (see below Eq. 29), as a
function of $m_{1}$, fixing $\delta=5\ \text{eV}$, $m_{A^{\prime}}/m_{1}=3$,
and $\alpha^{\prime}=0.5$. We have computed the $J$-factor by integrating over
the polar angles $0^{\circ}\leq\theta\leq 20^{\circ}$.
$J_{\text{iDM}}/J_{\text{NFW}}=1/4$ indicates that iDM coannihilations yield
the same flux as NFW-distributed Dirac DM annihilating with the same cross
section. For the red or blue lines, we assume an initial $\chi_{1,2}$
distribution described by an NFW profile or a $10\ \text{kpc}$ core Burkert
profile, respectively. The upper boundary of each shaded band labeled
“unvirialized” gives the $J$-factor assuming $\chi_{2}$ particles maintain the
same radius as their $\chi_{1}$ progenitors, whereas for the lower boundary of
each band labeled “revirialized” we assume that the $\chi_{2}$ population
eventually settles into the same type of Galactic distribution as $\chi_{1}$
with normalization fixed to conserve particle number. Also shown as dotted
lines are the contributions from the initial primordial density of excited
states, i.e., taking $\rho_{1,2}$ according to their initial values in Eq. 24.

Self-scattering processes can lead to significant modifications to the
$\chi_{1,2}$ density profiles. For instance, if the $\chi_{2}$ population is
sufficiently abundant, elastic scattering
$\chi_{1}\chi_{2}\to\chi_{1}\chi_{2}$ can result in a cored central density,
as typically found in models of self-interacting DM (SIDM) Tulin and Yu
(2018). However, the novel upscattering process
$\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$ alters the $\chi_{2}$ kinematics to
counteract such core formation, since it increases the likelihood for such
particles to cluster in the inner Galaxy. This is discussed briefly in
Appendix C and confirmed in the recent $N$-body simulation of Ref. O’Neil  _et
al._ (2022), which used model parameters comparable to the smallest DM masses
and mass-splittings that we consider in this work. Fully modeling all of these
effects requires a dedicated $N$-body simulation in our parameter space of
interest, which is beyond the scope of this work. Thus, in this section we
approximate the impact of these processes on the iDM $J$-factor by considering
a wide range of possible density profiles.

Figure 6: The effective annihilation cross section $\langle\sigma
v\rangle_{\text{eff}}$ (see Eq. 29), as a function of $m_{1}$ and for various
choices of $\delta$, fixing $m_{A^{\prime}}/m_{1}=3$, $\alpha^{\prime}=0.5$,
and $\epsilon$ such that the freeze-out abundance of $\chi_{1,2}$ agrees with
the observed DM density. The dotted lines show only the contribution from the
cosmological relic population of $\chi_{2}$ parameterized by $f_{2}$ in Eq.
13, whereas the solid lines also include the additional contribution from
Galactic upscattering in the Milky Way (assuming an initial NFW profile for
$\chi_{1}$ and $\chi_{2}$ and that the upscattered $\chi_{2}$ particles stay
fixed to the radii at which they are produced). Figure 7: Contours of the
total effective cross section $\langle\sigma v\rangle_{\text{eff}}$ (with
contributions from both the primordial and Galactic upscattered $\chi_{2}$
populations), as a function of $m_{1}$ and $\delta$, with the same model
assumptions as in Fig. 6. Also shown in gray are regions excluded by
observations of the CMB Aghanim  _et al._ (2020) and indirect detection
searches of the Milky Way Essig  _et al._ (2013); Boudaud  _et al._ (2017);
Cirelli  _et al._ (2021, 2023).

To determine $J_{\text{iDM}}$, we use $\rho_{1}$ and $\rho_{2}$ as calculated
in Sec. IV. Fig. 5 shows the ratio $J_{\text{iDM}}/J_{\text{NFW}}$ as a
function of $m_{1}$ for $\delta=5\ \text{eV}$. Here, we compute the $J$-factor
after integrating over a range of polar angles corresponding to
$0^{\circ}\leq\theta\leq 20^{\circ}$ and assuming either an initial NFW or
Burkert $\chi_{1}$ profile. In the latter case, the overall flux of the signal
is suppressed by our choice of the largest core size typically considered for
a Burkert profile, corresponding to $\sim 10\ \text{kpc}$ Nesti and Salucci
(2013) (this is meant to illustrate the maximal effect that coring from self-
scattering can have on $J_{\text{iDM}}$). Under each of these profile
assumptions, we investigate two further possibilities, corresponding to either
a $\chi_{2}$ population that tracks the same radius at which it is produced or
instead one that readjusts to the same type as the original $\chi_{1}$
profile, fixing the normalization to conserve total $\chi_{2}$ number (see the
discussion in Sec. IV). These two choices bracket the range of
$J_{\text{iDM}}$ values shown in Fig. 5. Also shown as dotted lines are the
contributions to the $J$-factor solely from the initial primordial density of
excited states, i.e., taking $\rho_{1,2}$ according to their initial values in
Eq. 24. Thus, we see that uncertainties regarding the effects of gravitational
dynamics can significantly impact the estimated flux of annihilation products,
at the level of a couple orders of magnitude within our parameter space of
interest.

In the remainder of this study, we therefore adopt the more optimistic
assumptions, corresponding to an initial NFW profile for $\chi_{1}$ and a
population of upscattered $\chi_{2}$ excited states that tracks the radial
position of their $\chi_{1}$ progenitors. However, we note that neither of the
bands shown in Fig. 5 take into account that the upscattered $\chi_{2}$
particles generically have less energy777Although the total energy is
conserved in these reactions, the $\chi_{2}$ produced from upscattering have
less combined kinetic and potential energy than their $\chi_{1}$ progenitors.
than their progenitors, which may lead to greater $\chi_{2}$ densities near
the Galactic Center and thereby enhance the $J$-factor relative to what is
shown here. We expand upon this phenomenon briefly in Appendix C and leave a
more detailed study to future work.

Figure 8:  Left: The effective cross section $\langle\sigma
v\rangle_{\text{eff}}$ (as in Eq. 29 and Fig. 6), as a function of $m_{1}$ for
various choices of $\delta$, fixing $m_{A^{\prime}}/m_{1}=3$,
$\alpha^{\prime}=0.5$, and $\epsilon$ such that the freeze-out abundance of
$\chi_{1,2}$ agrees with the observed DM density ($f_{\chi}=1$). The initial
distributions for $\chi_{1}$ and $\chi_{2}$ are taken to be NFW profiles, and
we assume that the upscattered $\chi_{2}$ stay fixed to the radii at which
they are produced. Three representative projections of proposed gamma-ray
satellites are shown as solid, dashed, and dotted black lines, including
e-ASTROGAM Tavani  _et al._ (2018), AMEGO Caputo  _et al._ (2022), and MAST
Dzhatdoev and Podlesnyi (2019), respectively. Also shown as shaded gray
regions are constraints from the CMB Aghanim  _et al._ (2020), indirect
detection searches of the Milky Way Essig  _et al._ (2013); Boudaud  _et al._
(2017); Cirelli  _et al._ (2021, 2023), and DM self-interactions in galaxy
clusters (labeled “SIDM”). The lower bound on $m_{1}$, as derived from the
successful predictions of BBN under the assumption of a standard cosmology
Giovanetti  _et al._ (2022), is shown as a vertical gray dotted line. The
approximate rate for $p$-wave annihilating thermal DM is shown as a dashed
pink line. Direct detection limits from iDM downscattering constrain
splittings larger than those considered here, $\delta\gtrsim 200\ \text{eV}$
Baryakhtar  _et al._ (2022); Carrillo González and Toro (2022). Right: Same as
the left panel, but instead fixing $\epsilon$ such that the freeze-out
abundance of $\chi_{1,2}$ constitutes only 10% of the total DM abundance
($f_{\chi}=0.1$). In this case, constraints on self-interacting DM do not
apply.

##  VI Results

In Fig. 6, we show $\langle\sigma v\rangle_{\text{eff}}$ (see Eq. 29) as a
function of $m_{1}$ and for various choices of $\delta$; the dotted lines show
the contribution from only the initial primordial population of excited
states. For sufficiently large splittings, upscattering in the Galaxy is
exponentially rare, such that the primordial population dominates the flux of
coannihilations at late times. From Eq. 20, for a fixed value of $\delta$,
upscattering is exponentially suppressed at low masses and power-law
suppressed at high masses, resulting in an optimal range of DM masses for the
Galactic contribution to $\langle\sigma v\rangle_{\text{eff}}$. On the other
hand, for the primordial population, $\langle\sigma v\rangle_{\text{eff}}$
grows with increasing mass due to the enhancement in $f_{2}$ (see the
discussion in Sec. III.4). This is evident in Fig. 6, which shows that for
splittings $\delta\lesssim 30\ \text{eV}$ and intermediate masses $m_{1}\sim
1\ \text{MeV}-100\ \text{MeV}$, coannihilations arising from the Galactic
upscattered $\chi_{2}$ particles dominantly contribute to $\langle\sigma
v\rangle_{\text{eff}}$. Note that for the primordial component, $f_{2}\simeq
1/2$ for $m_{1}\gtrsim 100\ \text{MeV}$, such that late-time annihilations are
unsuppressed and $\langle\sigma v\rangle_{\text{eff}}$ asymptotes to the
standard freeze-out value of $\langle\sigma v\rangle_{\text{eff}}\simeq
5\times 10^{-26}\ \text{cm}^{3}\,\text{s}^{-1}$ for sub-GeV masses Steigman
_et al._ (2012).

Fig. 7 displays $\langle\sigma v\rangle_{\text{eff}}$ in the two-dimensional
parameter space spanned by $m_{1}$ and $\delta$. In light gray we show regions
of parameter space excluded by observations of the CMB (as discussed in Sec.
III.4). Also shown in dark gray are constraints derived from existing
observations of the Milky Way. In particular, these consist of limits from
gamma-ray observations by INTEGRAL and COMPTEL Essig  _et al._ (2013), cosmic-
ray measurements by the Voyager 1 spacecraft Boudaud  _et al._ (2017), and
X-ray satellites Cirelli  _et al._ (2021, 2023). This region of parameter
space corresponds to scenarios in which $\chi_{2}$ efficiently reenters
equilibrium in the Galaxy, thus sourcing an annihilation rate that is
comparable to typical $s$-wave annihilations of thermal DM.

Proposed gamma-ray satellites will be able to significantly extend this
sensitivity to larger $\delta$, corresponding to smaller $\langle\sigma
v\rangle_{\text{eff}}$. This is evident in Fig. 8, which shows as colored
lines the predicted sizes of $\langle\sigma v\rangle_{\text{eff}}$ as a
function of mass, for various choices of $\delta$. In the left- and right-
panels, we fix $f_{\chi}=1$ and $f_{\chi}=0.1$, respectively. Also shown as a
dashed pink line is the approximate rate for $p$-wave thermal DM assuming
$f_{\chi}=1$. Three representative projections of proposed gamma-ray
satellites are shown as solid, dashed, and dotted black lines, including
e-ASTROGAM Tavani  _et al._ (2018), AMEGO Caputo  _et al._ (2022), and MAST
Dzhatdoev and Podlesnyi (2019), respectively. In particular, we have taken the
e-ASTROGAM projections from Ref. Bartels  _et al._ (2017) and the AMEGO and
MAST projections from Ref. Coogan  _et al._ (2021). We note that Ref. Coogan
_et al._ (2021) adopted an Einasto DM profile; we have rescaled their
projections appropriately to account for the smaller $J$-factor of an NFW
profile. Such instruments have the capability to probe a wide range of
currently unexplored thermal iDM parameter space, with annihilation rates
orders of magnitude larger than standard $p$-wave annihilating DM. Although
not shown in Fig. 8, the GECCO Orlando  _et al._ (2022) and COSI Zoglauer  _et
al._ (2021) instruments may also have competitive sensitivity in this
parameter space.

In addition to bounds from the CMB and existing indirect detection searches of
the Milky Way, for $f_{\chi}=1$ we also display limits from considerations of
self-interacting DM (SIDM). In particular, the large relative DM velocity
$v\sim 10^{-2}$ in observed galaxy clusters, such as the Bullet cluster
Robertson  _et al._ (2017), implies that DM as light as $m_{1}\sim 1\
\text{MeV}$ has enough kinetic energy to upscatter
$\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$ for splittings as large as $\delta\sim
100\ \text{eV}$. The region labeled “SIDM” in the left panel of Fig. 8
corresponds to parameter space where upscattering occurs at a rate of
$\sigma_{1\to 2}/m_{1}\gtrsim 10\ \text{cm}^{2}/\text{g}$. In this case, such
processes could lead to large deviations in the density profiles of merging
galaxy clusters. However, note that since typical studies of SIDM consider
elastic scattering, qualitative differences could arise for mass-splittings
comparable to the DM kinetic energy, such that properly recasting this bound
in the iDM parameter space necessitates a more careful treatment. Moreover,
since observations limit the total mass lost in such systems to be less than
$\sim 10\%$, iDM self-interactions are completely unconstrained if they make
up less than $\sim 10\%$ of the DM Markevitch  _et al._ (2004). This scenario
is considered in the right-panel of Fig. 8, where we instead consider a DM
subcomponent with $f_{\chi}=0.1$. Also note that in this case, the asymptotic
behavior of $\langle\sigma v\rangle_{\text{eff}}$ at high masses following
from Eqs. 10 and 29 is given by $\langle\sigma v\rangle_{\text{eff}}\propto
f_{\chi}^{2}\,\langle\sigma v\rangle_{\text{ann}}\propto f_{\chi}$.

For $m_{1}\lesssim 10\ \text{MeV}$, the thermal density of dark sector
particles can alter the expansion rate near temperatures of $T\sim 1\
\text{MeV}$, leading to predictions for the light element abundances that are
in conflict with standard Big Bang nucleosynthesis (BBN). In particular,
assuming a standard cosmological history, measurements of the CMB and inferred
primordial nuclear abundances exclude $m_{1}\lesssim 8\ \text{MeV}$ for the
dark sectors considered in this work Giovanetti  _et al._ (2022). This
constraint is shown as a dotted gray line in Fig. 8, but can be ameliorated if
additional new physics modifies the photon-neutrino temperature ratio or
introduces additional relativistic species at early times. We leave a more
complete model investigation along these lines to future work.

We conclude this section with a brief discussion on the sensitivity of direct
detection searches for DM scattering. In this work, we have refrained from
considering mass-splittings larger than $\delta=200\ \text{eV}$, since for
$200\ \text{eV}\lesssim\delta\lesssim 2m_{e}$, downscattering of primordial
excited states off of electrons, $\chi_{2}e\to\chi_{1}e$, is highly
constrained by existing XENON1T, SuperCDMS, and CRESST data Baryakhtar  _et
al._ (2022); Carrillo González and Toro (2022). Instead, for $\delta\lesssim
10^{-6}\,m_{1}$, upscattering in the form of $\chi_{1}e\to\chi_{2}e$ can be
searched for in low-threshold targets, but the limited exposure of current
experiments is insufficient to probe the thermal parameter space investigated
here Berlin  _et al._ (2019). Provided that low-backgrounds can be maintained
with larger exposures, future projections for the OSCURA Collaboration (2022)
and SuperCDMS Collaboration (2023) experiments could provide sensitivity to
such parameter space.

##  VII Discussion

We have explored a class of cosmologically-viable DM models in which the same
interaction that fixes the relic abundance also generates MeV-GeV gamma-rays
in the Galaxy that are detectable across a wide range of parameter space. In
this scenario, visible DM annihilations proceed efficiently in the early
Universe through the coannihilation between a lighter $\chi_{1}$ and heavier
$\chi_{2}$ state, but are exponentially suppressed near the time of
recombination as DM self-scattering continues to thermally deplete the
$\chi_{2}$ density. At much later times, the increased kinetic energy of DM
particles in the Milky Way enables $\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$
upscattering to efficiently regenerate a population of long-lived excited
states. Their resulting coannihilations give rise to low-energy gamma-rays
that are within reach of future telescopes, such as e-ASTROGAM Tavani  _et
al._ (2018), AMEGO Caputo  _et al._ (2022), and MAST Dzhatdoev and Podlesnyi
(2019). A novel feature of this scenario is that upscattering, followed by
coannihilation, is exponentially-dependent on the DM temperature, which peaks
at intermediate radii in the Milky Way. As a result, for certain DM mass-
splittings, there exists a radial shell within which $\chi_{1}$ and $\chi_{2}$
reenter chemical equilibrium over the age of the Galaxy.

Notably, such theories of inelastic DM serve as an important model-target for
accelerator and direct detection searches Krnjaic  _et al._ (2022). Hence,
indirect detection may play a key role in enabling model discrimination if
signals of sub-GeV DM arise (or fail to arise) elsewhere. For instance,
observing a signal in gamma-rays and in low-energy accelerator experiments,
such as LDMX Åkesson  _et al._ (2018); Berlin and Blinov (2018) or Belle-II
Campajola (2021), would provide evidence in favor of this scenario against
other models involving, e.g., $p$-wave annihilations (as is the case for
Majorana fermions or complex scalars annihilating to the SM through an
intermediate dark photon Berlin and Blinov (2018)).

In regards to direct detection searches, future experiments such as OSCURA
Collaboration (2022) and SuperCDMS Collaboration (2023) may ultimately be
sensitive to such models, either from downscattering $\chi_{2}e\to\chi_{1}e$
or upscattering $\chi_{1}e\to\chi_{2}e$, depending on the mass-
splitting.888Note that direct detection projections for inelastic (i.e.,
pseudo-Dirac) DM often assume that the leading signal arises from loop-induced
elastic scattering, which is beyond the reach of upcoming experiments.
However, as shown in Refs. Baryakhtar  _et al._ (2022); Carrillo González and
Toro (2022), downscattering can lead to larger signals, depending on the local
density of $\chi_{2}$.  Since in such models there is a unique prediction for
the relative strength of scattering signals and the Galactic annihilation
flux, mutually consistent signals from all three processes would constitute a
smoking gun discovery for this scenario and allow the determination of the
underlying model parameters.

In addition to what we have investigated in this work, there remain open
questions regarding cosmological and astrophysical predictions of these
models. Below, we highlight some future directions warranting further study
along these lines.

Galactic Dynamics: As noted above, our treatment of
$\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$ Galactic upscattering is based on an
analytical model of the Milky Way, which makes several idealized assumptions.
To quantify the theoretical uncertainty associated with our modeling, we have
considered a wide range of possibilities for the $\chi_{2}$ distribution in
the Galaxy (see Fig. 5). Nevertheless, there are several important effects
that cannot be easily estimated using our approach. For example, we separately
consider scenarios in which the upscattered $\chi_{2}$ population either
redistributes into the same halo profile as the $\chi_{1}$ population or,
alternatively, tracks the same radius where it was originally produced.
However, both possibilities ignore the fact that the $\chi_{2}$ are less
energetic than the $\chi_{1}$, due to upscattering. A brief discussion of this
point is given in Appendix C. A proper treatment of this effect is beyond the
scope of this paper and calls for a dedicated $N$-body numerical simulation,
along the lines of Ref. O’Neil  _et al._ (2022) but geared more towards the
specific parameter space investigated in our work. Since including additional
$\chi_{2}$ clustering would likely increase the Galactic $J$-factor, the
analytical projections in this paper are conservative.

Cosmological Evolution: Our calculation of the primordial $\chi_{2}$ fraction
adopts the instantaneous decoupling approximation in which $\chi_{2}$ number-
changing processes ($\chi_{1}e\to\chi_{2}e$ and
$\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$) maintain chemical equilibrium until
their rates fall below the Hubble scale. While this suffices for an order of
magnitude estimate of the cosmological $\chi_{2}$ number density, a more
accurate treatment requires numerically solving a system of Boltzmann
equations for $\chi_{1,2}$. We leave such an analysis to future work.

Cluster & Dwarf Comparison: For a given mass-splitting, a minimum amount of
kinetic energy is required to inelastically upscatter, and thus such processes
are highly sensitive to the effective DM temperature. As a result, the types
of signals discussed here may vary drastically across astrophysical systems,
such as dwarf galaxies and galaxy clusters, whose velocity dispersions differ
considerably from that of the Milky Way. In particular, for the parameter
space investigated in this work, such signals are exponentially-suppressed in
dwarf galaxies due to the smaller DM velocities involved McConnachie (2012);
O’Neil  _et al._ (2022). Conversely, it would be especially interesting to
study the sensitivity of proposed low-energy gamma-ray telescopes to signals
originating from galaxy clusters Lisanti  _et al._ (2018), as their sizeable
velocity dispersions Damsted  _et al._ (2023) could provide access to models
involving much larger mass-splittings. We defer a full treatment of these
issues to future work.

## Acknowledgements

We would like to thank Saniya Heeba, Dan Hooper, Katelin Schutz, and Oren
Slone for useful conversations. Fermilab is operated by the Fermi Research
Alliance, LLC under Contract DE-AC02-07CH11359 with the U.S. Department of
Energy. This material is based upon work supported by the U.S. Department of
Energy, Office of Science, National Quantum Information Science Research
Centers, Superconducting Quantum Materials and Systems Center (SQMS) under
contract number DE-AC02-07CH11359 and by the Kavli Institute for Cosmological
Physics at the University of Chicago through an endowment from the Kavli
Foundation and its founder Fred Kavli. This work was completed in part at the
Perimeter Institute. Research at Perimeter Institute is supported in part by
the Government of Canada through the Department of Innovation, Science and
Economic Development Canada and by the Province of Ontario through the
Ministry of Colleges and Universities.

##  Appendix A $\text{DM}-\text{SM}$ Kinetic Decoupling

Here, we provide discussion regarding $\text{DM}-\text{SM}$ scattering and
kinetic decoupling. In the limit that $m_{A^{\prime}}\gg\sqrt{s}\gg\delta$ and
$m_{1}\gg m_{e},T$, the differential cross section for $\chi_{2}e\to\chi_{1}e$
is given by

$\frac{d\sigma_{\chi
e}}{dt}\simeq\frac{\pi\alpha\alpha^{\prime}\epsilon^{2}}{2m_{A^{\prime}}^{4}m_{1}^{2}\,p_{e}^{2}}\,\Big{[}2(s-m_{1}^{2}-m_{e}^{2})^{2}+2st+t^{2}\Big{]}~{},$
(30)

where $s$ and $t$ are the usual Mandelstam variables and $p_{e}$ is the
electron momentum. Following Refs. Gondolo  _et al._ (2012); Bertoni  _et al._
(2015), the rate for maintaining kinetic equilibrium is then given by

$\Gamma_{\text{kd}}\simeq-\frac{1}{3m_{1}T}\,\int\frac{d^{3}p_{e}}{(2\pi)^{3}}~{}f_{e}(1-f_{e})\,v_{e}\,\int_{-4p_{e}^{2}}^{0}dt~{}t\,\frac{d\sigma_{\chi
e}}{dt}~{},$ (31)

where $f_{e}=(e^{E_{e}/T}+1)^{-1}$ is the phase-space distribution of
electrons with energy $E_{e}$, and $v_{e}$ is the electron velocity. In the
limit $T\gg m_{e}$ or $T\ll m_{e}$, the above expression reduces to Eq. 11.
Note that schematically, the above rate is $\Gamma_{\text{kd}}\sim
n_{e}\sigma_{\chi e}v_{e}/N_{\text{scatt}}$, where $N_{\text{scatt}}\sim
m_{1}T/q^{2}\gg 1$ is the number of scatters, each exchanging momentum $q$,
needed to transfer a total amount comparable to the DM momentum
$\sim\sqrt{m_{1}T}$.

##  Appendix B $\text{DM}-\text{DM}$ Chemical Decoupling

Here, we provide more discussion regarding DM downscattering
$\chi_{2}\chi_{2}\to\chi_{1}\chi_{1}$. To begin, we calculate the cross
section in the limit that $m_{A^{\prime}}\gg\sqrt{s}$ to be

$\sigma_{2\to
1}\simeq\frac{2\pi\alpha^{\prime\,2}}{3sm_{A^{\prime}}^{4}}\,\sqrt{\frac{s-4m_{1}^{2}}{s-4m_{2}^{2}}}~{}\Big{[}7s^{2}-4s\,(4m_{1}^{2}+4m_{2}^{2}+3m_{1}m_{2})+76m_{1}^{2}m_{2}^{2}\Big{]}~{}.$
(32)

To evaluate the thermal average $\langle\sigma v\rangle_{2\to 1}$, we use Eq.
3.8 of Ref. Gondolo and Gelmini (1991). In the limit that $\delta\ll
T_{\chi}\ll m_{1}$ or $T_{\chi}\ll\delta\ll m_{1}$, we arrive at Eq. 21. The
corresponding downscattering rate of Eq. 12 is then given by
$\Gamma_{\chi}=n_{2}\,\langle\sigma v\rangle_{2\to 1}$, where
$n_{2}=e^{-\delta/T_{\chi}}\,n_{1}$.

##  Appendix C Upscattering in the Milky Way Halo

Since upscattering $\chi_{1}\chi_{1}\to\chi_{2}\chi_{2}$ converts kinetic
energy into mass energy, it causes particles to fall more deeply into
gravitational wells, enhancing the central density of halos, analogous to the
dynamics associated with dissipative reactions. This was noted in the detailed
numerical simulations of Ref. O’Neil  _et al._ (2022), where it was found that
enhanced central densities from upscattering can counteract coring from
general scattering-induced heat transport throughout the halo. Here, we
provide a brief semi-analytic discussion of the central density enhancement
from upscattering, leaving a more dedicated analysis to future work.

We denote the total energy of an initial $\chi_{1}$ particle as $E=m_{1}+T+U$,
where $T$ and $U$ are its kinetic and gravitational potential energy,
respectively. For convenience, let us rewrite this as $E=m_{1}+\tilde{E}$,
where $\tilde{E}\equiv T+U$. Using the virial theorem, $T=-\frac{1}{2}U$, we
can rewrite $\tilde{E}$ solely in terms of the gravitational energy at radius
$r$,

$\tilde{E}=\frac{1}{2}\,U=-\frac{GM_{\text{enc}}(r)\,m_{1}}{2r}~{},$ (33)

where $M_{\text{enc}}(r)$ is the total Galactic mass enclosed within $r$.
Using that the total energy $E$ is conserved as $\chi_{1}$ upscatters and
converts into the slightly heavier state $\chi_{2}$, we can relate the change
in $\tilde{E}$ to the change in mass $m_{\chi}$ of the particle,

$\frac{d\tilde{E}}{dt}=-\frac{dm_{\chi}}{dt}=-\Gamma_{1\to 2}~{}\delta~{},$
(34)

where $\Gamma_{1\to 2}=n_{1}\,\langle\sigma v\rangle_{1\to 2}$ is the
upscattering rate per $\chi_{1}$ particle. Alternatively, the change in
$\tilde{E}$ can be rewritten using the chain rule as

$\frac{d\tilde{E}}{dt}=\frac{d\tilde{E}}{dr}\,\frac{dr}{dt}=\frac{GM_{\text{enc}}(r)\,m_{1}}{2r^{2}}\,\frac{dr}{dt}~{},$
(35)

where we used Eq. 33 in the last equality. Equating Eqs. 34 and 35, we have
that the change in radius from upscattering is given by

$\frac{dr}{dt}\simeq-2\,\Gamma_{1\to
2}~{}\frac{\delta}{m_{1}}\,\frac{r^{2}}{GM_{\text{enc}}(r)}~{}.$ (36)

Given a general profile for the enclosed mass $M_{\text{enc}}(r)$ and an
initial radius $r_{i}=r(0)$ at $t=0$, Eq. 36 can be solved for the radius
today $r(t_{\text{mw}})$. This can then be inverted to obtain $r_{i}(r)$.
Since upscattering does not change total particle number, we then relate the
initial number density profile at early times $n_{i}(r_{i})$ to the profile
today $n(r)$ via $n_{i}(r_{i})\,r_{i}^{2}=n(r)\,r^{2}$. Thus, to leading order
in the mass-splitting, the mass density profiles are similarly related
$\rho_{i}(r_{i})\,r_{i}^{2}\simeq\rho(r)\,r^{2}$. The final form of the mass
density profile is then determined by

$\rho(r)\simeq\bigg{(}\frac{r_{i}(r)}{r}\bigg{)}^{2}\,\rho_{i}\big{[}r_{i}(r)\big{]}~{},$
(37)

where the function $r_{i}(r)$ is determined from Eq. 36, as described above.
We have found that Eq. 37 accurately captures the salient features found in
the “endothermic up-scattering” simulation of Ref. O’Neil  _et al._ (2022).
However, note that the above analysis does not incorporate the full set of
dynamics, as it only allows for $\chi_{1}\to\chi_{2}$ conversion and not the
reverse process.

## References

  * Krnjaic  _et al._ (2022) G. Krnjaic _et al._ , “A Snowmass Whitepaper: Dark Matter Production at Intensity-Frontier Experiments,” (2022), arXiv:2207.00597 [hep-ph] .
  * Essig  _et al._ (2022) Rouven Essig _et al._ , “Snowmass2021 Cosmic Frontier: The landscape of low-threshold dark matter direct detection in the next decade,” in  _Snowmass 2021_ (2022) arXiv:2203.08297 [hep-ph] .
  * Berlin  _et al._ (2019) Asher Berlin, Nikita Blinov, Gordan Krnjaic, Philip Schuster, and Natalia Toro, “Dark Matter, Millicharges, Axion and Scalar Particles, Gauge Bosons, and Other New Physics with LDMX,” Phys. Rev. D 99, 075001 (2019), arXiv:1807.01730 [hep-ph] .
  * Padmanabhan and Finkbeiner (2005) Nikhil Padmanabhan and Douglas P. Finkbeiner, “Detecting dark matter annihilation with CMB polarization: Signatures and experimental prospects,” Phys. Rev. D 72, 023508 (2005), arXiv:astro-ph/0503486 .
  * Slatyer (2016) Tracy R. Slatyer, “Indirect dark matter signatures in the cosmic dark ages. I. Generalizing the bound on s-wave dark matter annihilation from Planck results,” Phys. Rev. D 93, 023527 (2016), arXiv:1506.03811 [hep-ph] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Feng and Smolinsky (2017) Jonathan L. Feng and Jordan Smolinsky, “Impact of a resonance on thermal targets for invisible dark photon searches,” Phys. Rev. D 96, 095022 (2017), arXiv:1707.03835 [hep-ph] .
  * Bernreuther  _et al._ (2021) Elias Bernreuther, Saniya Heeba, and Felix Kahlhoefer, “Resonant sub-GeV Dirac dark matter,” JCAP 03, 040 (2021), arXiv:2010.14522 [hep-ph] .
  * Brahma  _et al._ (2023) Nirmalya Brahma, Saniya Heeba, and Katelin Schutz, “Resonant Pseudo-Dirac Dark Matter as a Sub-GeV Thermal Target,” (2023), arXiv:2308.01960 [hep-ph] .
  * Griest and Seckel (1991) Kim Griest and David Seckel, “Three exceptions in the calculation of relic abundances,” Phys. Rev. D 43, 3191–3203 (1991).
  * D’Agnolo and Ruderman (2015) Raffaele Tito D’Agnolo and Joshua T. Ruderman, “Light Dark Matter from Forbidden Channels,” Phys. Rev. Lett. 115, 061301 (2015), arXiv:1505.07107 [hep-ph] .
  * D’Agnolo  _et al._ (2021) Raffaele Tito D’Agnolo, Di Liu, Joshua T. Ruderman, and Po-Jen Wang, “Forbidden dark matter annihilations into Standard Model particles,” JHEP 06, 103 (2021), arXiv:2012.11766 [hep-ph] .
  * Bartels  _et al._ (2017) Richard Bartels, Daniele Gaggero, and Christoph Weniger, “Prospects for indirect dark matter searches with MeV photons,” JCAP 05, 001 (2017), arXiv:1703.02546 [astro-ph.HE] .
  * Coogan  _et al._ (2021) Adam Coogan, Logan Morrison, and Stefano Profumo, “Precision gamma-ray constraints for sub-GeV dark matter models,” JCAP 08, 044 (2021), arXiv:2104.06168 [hep-ph] .
  * Tucker-Smith and Weiner (2001) David Tucker-Smith and Neal Weiner, “Inelastic dark matter,” Phys. Rev. D 64, 043502 (2001), arXiv:hep-ph/0101138 .
  * Arkani-Hamed  _et al._ (2009) Nima Arkani-Hamed, Douglas P. Finkbeiner, Tracy R. Slatyer, and Neal Weiner, “A Theory of Dark Matter,” Phys. Rev. D 79, 015014 (2009), arXiv:0810.0713 [hep-ph] .
  * Arkani-Hamed and Weiner (2008) Nima Arkani-Hamed and Neal Weiner, “LHC Signals for a SuperUnified Theory of Dark Matter,” JHEP 12, 104 (2008), arXiv:0810.0714 [hep-ph] .
  * Finkbeiner and Weiner (2007) Douglas P. Finkbeiner and Neal Weiner, “Exciting Dark Matter and the INTEGRAL/SPI 511 keV signal,” Phys. Rev. D 76, 083519 (2007), arXiv:astro-ph/0702587 .
  * Chen  _et al._ (2010) Fang Chen, James M. Cline, Anthony Fradette, Andrew R. Frey, and Charles Rabideau, “Exciting dark matter in the galactic center,” Phys. Rev. D 81, 043523 (2010), arXiv:0911.2222 [hep-ph] .
  * Finkbeiner and Weiner (2016) Douglas P. Finkbeiner and Neal Weiner, “X-ray line from exciting dark matter,” Phys. Rev. D 94, 083002 (2016), arXiv:1402.6671 [hep-ph] .
  * D’Eramo and Profumo (2018) Francesco D’Eramo and Stefano Profumo, “Sub-GeV Dark Matter Shining at Future MeV $\gamma$ -Ray Telescopes,” Phys. Rev. Lett. 121, 071101 (2018), arXiv:1806.04745 [hep-ph] .
  * Buckley and Profumo (2012) Matthew R. Buckley and Stefano Profumo, “Regenerating a Symmetry in Asymmetric Dark Matter,” Phys. Rev. Lett. 108, 011301 (2012), arXiv:1109.2164 [hep-ph] .
  * Holdom (1986) Bob Holdom, “Two U(1)’s and Epsilon Charge Shifts,” Phys. Lett. B 166, 196–198 (1986).
  * Baryakhtar  _et al._ (2022) Masha Baryakhtar, Asher Berlin, Hongwan Liu, and Neal Weiner, “Electromagnetic signals of inelastic dark matter scattering,” JHEP 06, 047 (2022), arXiv:2006.13918 [hep-ph] .
  * Carrillo González and Toro (2022) Mariana Carrillo González and Natalia Toro, “Cosmology and signals of light pseudo-Dirac dark matter,” JHEP 04, 060 (2022), arXiv:2108.13422 [hep-ph] .
  * Evans  _et al._ (2018) Jared A. Evans, Stefania Gori, and Jessie Shelton, “Looking for the WIMP Next Door,” JHEP 02, 100 (2018), arXiv:1712.03974 [hep-ph] .
  * Berlin  _et al._ (2014) Asher Berlin, Dan Hooper, and Samuel D. McDermott, “Simplified Dark Matter Models for the Galactic Center Gamma-Ray Excess,” Phys. Rev. D 89, 115022 (2014), arXiv:1404.0022 [hep-ph] .
  * Izaguirre  _et al._ (2016) Eder Izaguirre, Gordan Krnjaic, and Brian Shuve, “Discovering Inelastic Thermal-Relic Dark Matter at Colliders,” Phys. Rev. D 93, 063523 (2016), arXiv:1508.03050 [hep-ph] .
  * Berlin and Kling (2019) Asher Berlin and Felix Kling, “Inelastic Dark Matter at the LHC Lifetime Frontier: ATLAS, CMS, LHCb, CODEX-b, FASER, and MATHUSLA,” Phys. Rev. D 99, 015021 (2019), arXiv:1810.01879 [hep-ph] .
  * Berlin  _et al._ (2020) Asher Berlin, Patrick deNiverville, Adam Ritz, Philip Schuster, and Natalia Toro, “Sub-GeV dark matter production at fixed-target experiments,” Phys. Rev. D 102, 095011 (2020), arXiv:2003.03379 [hep-ph] .
  * Batell  _et al._ (2009) Brian Batell, Maxim Pospelov, and Adam Ritz, “Direct Detection of Multi-component Secluded WIMPs,” Phys. Rev. D 79, 115019 (2009), arXiv:0903.3396 [hep-ph] .
  * McMillan (2011) Paul J. McMillan, “Mass models of the Milky Way,” Mon. Not. Roy. Astron. Soc. 414, 2446–2457 (2011), arXiv:1102.4340 [astro-ph.GA] .
  * Kaplinghat  _et al._ (2014) Manoj Kaplinghat, Ryan E. Keeley, Tim Linden, and Hai-Bo Yu, “Tying Dark Matter to Baryons with Self-interactions,” Phys. Rev. Lett. 113, 021302 (2014), arXiv:1311.6524 [astro-ph.CO] .
  * Navarro  _et al._ (1997) Julio F. Navarro, Carlos S. Frenk, and Simon D. M. White, “A Universal density profile from hierarchical clustering,” Astrophys. J. 490, 493–508 (1997), arXiv:astro-ph/9611107 .
  * Tulin and Yu (2018) Sean Tulin and Hai-Bo Yu, “Dark Matter Self-interactions and Small Scale Structure,” Phys. Rept. 730, 1–57 (2018), arXiv:1705.02358 [hep-ph] .
  * O’Neil  _et al._ (2022) Stephanie O’Neil _et al._ , “Endothermic self-interacting dark matter in Milky Way-like dark matter haloes,”  (2022), 10.1093/mnras/stad1850, arXiv:2210.16328 [astro-ph.GA] .
  * Essig  _et al._ (2013) Rouven Essig, Eric Kuflik, Samuel D. McDermott, Tomer Volansky, and Kathryn M. Zurek, “Constraining Light Dark Matter with Diffuse X-Ray and Gamma-Ray Observations,” JHEP 11, 193 (2013), arXiv:1309.4091 [hep-ph] .
  * Boudaud  _et al._ (2017) Mathieu Boudaud, Julien Lavalle, and Pierre Salati, “Novel cosmic-ray electron and positron constraints on MeV dark matter particles,” Phys. Rev. Lett. 119, 021103 (2017), arXiv:1612.07698 [astro-ph.HE] .
  * Cirelli  _et al._ (2021) Marco Cirelli, Nicolao Fornengo, Bradley J. Kavanagh, and Elena Pinetti, “Integral X-ray constraints on sub-GeV Dark Matter,” Phys. Rev. D 103, 063022 (2021), arXiv:2007.11493 [hep-ph] .
  * Cirelli  _et al._ (2023) Marco Cirelli, Nicolao Fornengo, Jordan Koechler, Elena Pinetti, and Brandon M. Roach, “Putting all the X in one basket: Updated X-ray constraints on sub-GeV Dark Matter,” JCAP 07, 026 (2023), arXiv:2303.08854 [hep-ph] .
  * Nesti and Salucci (2013) Fabrizio Nesti and Paolo Salucci, “The Dark Matter halo of the Milky Way, AD 2013,” JCAP 07, 016 (2013), arXiv:1304.5127 [astro-ph.GA] .
  * Tavani  _et al._ (2018) M. Tavani _et al._ (e-ASTROGAM), “Science with e-ASTROGAM: A space mission for MeV–GeV gamma-ray astrophysics,” JHEAp 19, 1–106 (2018), arXiv:1711.01265 [astro-ph.HE] .
  * Caputo  _et al._ (2022) Regina Caputo _et al._ , “All-sky Medium Energy Gamma-ray Observatory eXplorer mission concept,” J. Astron. Telesc. Instrum. Syst. 8, 044003 (2022), arXiv:2208.04990 [astro-ph.IM] .
  * Dzhatdoev and Podlesnyi (2019) Timur Dzhatdoev and Egor Podlesnyi, “Massive Argon Space Telescope (MAST): A concept of heavy time projection chamber for $\gamma$-ray astronomy in the 100 MeV–1 TeV energy range,” Astropart. Phys. 112, 1–7 (2019), arXiv:1902.01491 [astro-ph.HE] .
  * Giovanetti  _et al._ (2022) Cara Giovanetti, Mariangela Lisanti, Hongwan Liu, and Joshua T. Ruderman, “Joint Cosmic Microwave Background and Big Bang Nucleosynthesis Constraints on Light Dark Sectors with Dark Radiation,” Phys. Rev. Lett. 129, 021302 (2022), arXiv:2109.03246 [hep-ph] .
  * Steigman  _et al._ (2012) Gary Steigman, Basudeb Dasgupta, and John F. Beacom, “Precise Relic WIMP Abundance and its Impact on Searches for Dark Matter Annihilation,” Phys. Rev. D 86, 023506 (2012), arXiv:1204.3622 [hep-ph] .
  * Orlando  _et al._ (2022) Elena Orlando _et al._ , “Exploring the MeV sky with a combined coded mask and Compton telescope: the Galactic Explorer with a Coded aperture mask Compton telescope (GECCO),” JCAP 07, 036 (2022), arXiv:2112.07190 [astro-ph.HE] .
  * Zoglauer  _et al._ (2021) Andreas Zoglauer _et al._ , “COSI: From Calibrations and Observations to All-sky Images,” (2021), arXiv:2102.13158 [astro-ph.IM] .
  * Robertson  _et al._ (2017) Andrew Robertson, Richard Massey, and Vincent Eke, “What does the Bullet Cluster tell us about self-interacting dark matter?” Mon. Not. Roy. Astron. Soc. 465, 569–587 (2017), arXiv:1605.04307 [astro-ph.CO] .
  * Markevitch  _et al._ (2004) Maxim Markevitch, A. H. Gonzalez, D. Clowe, A. Vikhlinin, L. David, W. Forman, C. Jones, S. Murray, and W. Tucker, “Direct constraints on the dark matter self-interaction cross-section from the merging galaxy cluster 1E0657-56,” Astrophys. J. 606, 819–824 (2004), arXiv:astro-ph/0309303 .
  * Collaboration (2022) Oscura Collaboration, “The oscura experiment,” (2022), arXiv:2202.10518 [astro-ph.IM] .
  * Collaboration (2023) SuperCDMS Collaboration, “A strategy for low-mass dark matter searches with cryogenic detectors in the supercdms snolab facility,” (2023), arXiv:2203.08463 [physics.ins-det] .
  * Åkesson  _et al._ (2018) Torsten Åkesson _et al._ (LDMX), “Light Dark Matter eXperiment (LDMX),” (2018), arXiv:1808.05219 [hep-ex] .
  * Berlin and Blinov (2018) Asher Berlin and Nikita Blinov, “Thermal Dark Matter Below an MeV,” Phys. Rev. Lett. 120, 021801 (2018), arXiv:1706.07046 [hep-ph] .
  * Campajola (2021) Marcello Campajola (Belle-II), “Dark Sector first results at Belle II,” Phys. Scripta 96, 084005 (2021), arXiv:2208.01101 [hep-ex] .
  * McConnachie (2012) Alan W. McConnachie, “The observed properties of dwarf galaxies in and around the Local Group,” Astron. J. 144, 4 (2012), arXiv:1204.1562 [astro-ph.CO] .
  * Lisanti  _et al._ (2018) Mariangela Lisanti, Siddharth Mishra-Sharma, Nicholas L. Rodd, and Benjamin R. Safdi, “Search for Dark Matter Annihilation in Galaxy Groups,” Phys. Rev. Lett. 120, 101101 (2018), arXiv:1708.09385 [astro-ph.CO] .
  * Damsted  _et al._ (2023) S. Damsted, A. Finoguenov, N. Clerc, I. Davalgaite, C. C. Kirkpatrick, G. A. Mamon, J. Ider Chitham, K. Kiiveri, J. Comparat, and C. Collins, “CODEX: Role of velocity substructure in the scaling relations of galaxy clusters,” (2023), arXiv:2307.08749 [astro-ph.CO] .
  * Gondolo  _et al._ (2012) Paolo Gondolo, Junji Hisano, and Kenji Kadota, “The Effect of quark interactions on dark matter kinetic decoupling and the mass of the smallest dark halos,” Phys. Rev. D 86, 083523 (2012), arXiv:1205.1914 [hep-ph] .
  * Bertoni  _et al._ (2015) Bridget Bertoni, Seyda Ipek, David McKeen, and Ann E. Nelson, “Constraints and consequences of reducing small scale structure via large dark matter-neutrino interactions,” JHEP 04, 170 (2015), arXiv:1412.3113 [hep-ph] .
  * Gondolo and Gelmini (1991) Paolo Gondolo and Graciela Gelmini, “Cosmic abundances of stable particles: Improved analysis,” Nucl. Phys. B 360, 145–179 (1991).