(see also Fig. 5). Firstly, the number of bound states at a given mass is
shown on the left of Fig. 6. We have argued that production of a bound state
with a given $S$ and $L>0$ will generically lead to production (via one or
more decays) of an $L=0$ bound state with the same $S$, so it unnecessary to
compute the rate to produce SM particles directly from those $L>0$ states. We
confirm this numerically on the right of Fig. 6: however measured, the
branching fraction is always greater than 90%. For instance, at 100 TeV, the
deepest $Q=0$ $L+S$ odd $p$-wave bound state decays to the deepest $s$-wave
state with 99.6% probability and to the second deepest $s$-wave state with a
likelihood of 0.4%. (The next most dominant transition is to the third deepest
$s$-wave bound state, and only occurs with a probability of 0.006%.)

We also discussed decays of the $S=1,\,L=0$ bound states, noting that they
generate only a power suppressed contribution to the photon endpoint spectrum.
Nevertheless, the contribution from the $S=1$ states can be comparable to that
from the $S=0$ states, after accounting for both their (enhanced) production
rate and (suppressed) endpoint spectrum from annihilation. However, in
practice this means that the contributions to the endpoint spectrum from both
the $S=0$ and $S=1$ bound states are suppressed compared to direct
annihilation (either due to power suppression or because of the small
formation rate).

For the aforementioned reasons, we focus on the $L=S=0$ bound states to
estimate the size of the (generally subdominant) bound-state contribution.
Doing so implies we can reuse our SCET results from the direct annihilation
case, with appropriate modification of the wavefunction factors. This choice
does mean that the total bound-state contribution to the endpoint photon
spectrum could increase by an $\mathcal{O}(1)$ factor in an improved
calculation (once the $S=1$ states are included). Because the bound-state
contribution to the endpoint spectrum is suppressed, this typically
corresponds to a percent-level theoretical uncertainty in the overall endpoint
spectrum, with larger uncertainties at specific mass points where the bound
state formation cross section is enhanced relative to direct annihilation. The
contribution to continuum photons (not near the endpoint) from bound state
formation can be markedly larger, compared to the effect on the endpoint
spectrum, as there is no power suppression in the continuum contribution from
annihilation of the $S=1$ bound states. We will discuss each of these
contributions to the total spectrum in the next section.

##  5 The Combined Photon Spectrum and Numerical Results

At this stage we have all ingredients required to determine the quintuplet
annihilation spectrum, including both direct annihilation and the contribution
of bound states. In this section we collect our results to determine the full
energy distribution of photons the quintuplet generates at the thermal mass of
13.6 TeV, but also for a wider range of masses. We will estimate the impact of
several uncertainties on our results, such as the residual theoretical
uncertainty on the NLL computations, but also from astrophysical uncertainties
such as the distribution of $v$ values and on the DM density in the inner
Galaxy. Finally, we will put these results together to estimate the
sensitivity of existing and upcoming IACTs to quintuplet DM.

###  5.1 Predictions for the spectrum and rate of photon production

A central goal of this work is to accurately determine the distribution of
photons that emerge when two SU(2) quintuplets annihilate. This spectrum forms
the signal template for telescopes searching for high energy photons, and
therefore is a central theoretical input. To achieve this, throughout we have
computed differential cross sections $d\sigma/dz$, both for the direct
annihilation in Eq. (31), and also for the bound state contribution by
combining the results of the previous sections with Eq. (1). For indirect
detection, observables are sensitive to $d\langle\sigma v\rangle/dz$. To begin
with we will assume the DM states are incident with a fixed $v=10^{-3}$,
revisiting the validity of this approximation in the next subsection. In order
to extract the shape of the photon distribution from the differential cross
section, it is common in indirect detection to introduce a photon spectrum
$dN/dE$, and our convention for doing so is the following,151515Further
discussion of the connection between spectra used in indirect detection and
the corresponding field theoretic quantities can be found in, for instance,
Ref. Bauer:2020jay .

$\displaystyle\frac{d\langle\sigma v\rangle}{dE}=\langle\sigma v\rangle_{\rm
line}\times\frac{dN}{dE}.$ (84)

This choice follows Ref. Baumgart:2017nsr , and implies that our spectrum is
normalized with respect to the _line cross section_ , $\langle\sigma
v\rangle_{\rm line}\equiv\langle\sigma
v\rangle_{\gamma\gamma}+\tfrac{1}{2}\langle\sigma v\rangle_{\gamma Z}$, which
is defined as the rate to produce two photons at exactly $E=M_{\chi}$. By
construction, $dN/dE$ will contain a contribution of exactly
$2\delta(E-M_{\chi})$ for the line, but it will also contain contributions
from endpoint photons, bound state decays, and continuum photons arising
primarily from the unstable particles the direct annihilations or bound state
decays can produce. For each of these latter components, however, their
additions to $dN/dE$ will be weighted by a branching fraction $\langle\sigma
v\rangle_{i}/\langle\sigma v\rangle_{\rm line}$, with $\langle\sigma
v\rangle_{i}$ the cross section for that particular contribution. The
rationale for anchoring our calculations to the line cross section is that
$\chi\chi\to\gamma\gamma$, which has a spectrum of exactly
$dN/dE=2\delta(E-M_{\chi})$, is a common experimental target, and therefore
there are a wide number of existing constraints on $\langle\sigma
v\rangle_{\rm line}$ which we can then directly compare with. Further
discussion of this point can be found in Ref. Baumgart:2017nsr .

Figure 7: The cross section for line photons, breaking down the contributions
from the direct annihilation and bound states. While generically the direct
annihilation dominates, for isolated masses near Sommerfeld peaks the bound
state contribution can be the leading one. For all masses lower than those
shown bound states are strictly subdominant.

The spectra of line and endpoint photons produced by decay of the bound states
is computed using the methods of Sec. 4.161616In this work we do not track the
much lower-energy photons radiated in bound state formation and transitions;
these signals have been studied for the quintuplet in Refs. Mitridate:2017izz
; Mahbubani:2020knq . For our bound state formation and transition
calculations, we include only states with $L=0,1,2$. Capture into $L=3$ and
higher states requires at least $L=2$ for the initial state, and we expect the
contributions from components of the initial state with high $L$ to be
suppressed at low velocities, by a factor that is parametrically
$(M_{\chi}v/m_{\scriptscriptstyle W})^{2L}$ (although for sufficiently high
masses, $M_{\chi}\gtrsim 100$ TeV, this suppression is lifted). It is also
worth noting that for essentially all the parameter space most relevant for
experimental searches with H.E.S.S, at $M_{\chi}<20$ TeV, we find that no
$L=3+$ bound states exist in the spectrum. We independently expect capture
into states with high principal quantum number $n$ (which is required for high
$L$) to be suppressed, as (1) the finite range of the potential means only a
limited number of states are bound at all, so unlike in unbroken gauge
theories there is no infinite tower of high-$n$ states, (2) capture into
weakly-bound states is suppressed by a phase space factor, and (3) analytic
approximations (App. F) suggest that we can expect the leading contribution to
the capture cross section to be exponentially suppressed for large $n$. In
practice, our numerical calculation expresses the bound states as a linear
combination of 30 basis states for each combination of $L$, $Q$ and
$(-1)^{L+S}$, allowing us to access up to 30 distinct bound states indexed by
different values of $n$ (although as we approach this upper bound we expect
the spectrum to become less accurate), and we include all these states in our
calculation. We have checked at sample mass points that our binding energies
and cross sections for capture into lower-$n$ states are not significantly
affected by doubling the number of basis states. For the reasons given above,
we generally expect the error due to the omission of higher-$n$ states to be
small.

Before showing the full distributions of line and endpoint photons, we can
already consider one measure of the importance of bound states to the
resulting photon signal: their contribution to $\langle\sigma
v\rangle_{\text{line}}$. This is shown in Fig. 7, where we separate the
contribution to the line from the direct annihilation to that of processes
involving an intermediate bound state. At this stage we only include bound-
state contributions that produce line photons, with energy essentially at
$M_{\chi}$. The figure makes clear a point already estimated earlier: direct
annihilation generally dominates the production of line photons at $E\simeq
M_{\chi}$ by $1-2$ orders of magnitude. However, the bound-state contribution
can be significant and even dominate at isolated mass points, for instance as
at $M_{\chi}=68.1~{}{\rm TeV}$, and therefore a reliable prediction at
arbitrary masses must include this contribution.

Figure 8: The quintuplet annihilation spectrum for two masses, the thermal
mass of 13.6 TeV (left), and a mass where bound state contributions are
appreciable, 68.1 TeV (right). For each mass we show results convolved with
the H.E.S.S. energy resolution (top) and unsmeared (below). The full spectrum
is broken into five individual components: the line, endpoint, bound state
line and endpoint, the direct annihilation continuum, and bound state
contribution to the continuum. Details of each are provided in the text.

Moving beyond the line, in Fig. 8 we show the full spectrum, broken down by
various contributions, for two masses: the thermal mass of
$M_{\chi}=13.6~{}{\rm TeV}$, and a mass where bound state contributions are
significant, $68.1~{}{\rm TeV}$. For each mass we show two versions of the
spectrum. In the lower panels, we show the unsmeared spectra, which is the
distribution of photons that emerge from the annihilations. (Note in this case
the line contribution is simply $dN/dE=2\delta(E-M_{\chi})$ and so represented
by a vertical line.) In the upper panels, we have convolved the raw spectra
with a finite experimental energy resolution in order to model what would
actually be seen at a realistic instrument. For this we take the energy
resolution of the H.E.S.S. telescope, determined from Ref. HESS:2013rld . In
detail, we fix the relative width $\Delta E/E$ to 0.17 and 0.11 for
$E=500~{}{\rm GeV}$ and $E=10~{}{\rm TeV}$, respectively, and then vary
logarithmically between these endpoints, freezing the ratio either side. From
this we compute $(dN/dE)_{\rm smeared}$ as $dN/dE$ convolved with a Gaussian
of width equal to the energy resolution.

In terms of these two notions of the spectra, Fig. 8 shows five contributions
to the photons distributions for the two masses. The first three of these are:
1. the direct annihilation line; 2. direct annihilation endpoint; and 3. the
bound state contribution to the line and endpoint. Again we see clearly a
point noted for the wino in Refs. Baumgart:2017nsr ; Baumgart:2018yed : the
endpoint contribution makes a considerable modification to the observed number
of photons with $E\sim M_{\chi}$, with the peak smeared spectra enhanced by
1.9 and 3.1 for $M_{\chi}=13.6$ and 68.1 TeV. The bound state contribution is
more modest: it is effectively negligible at the thermal mass, and a factor of
1.7 enhancement at 68.1 TeV, which again is a mass with an anomalously large
bound state contribution to the hard photon spectrum.

Beyond these three, we also show the contribution of two continuum sources,
both of which can generate lower energy photons. The first of these is the
continuum emission arising from direct annihilation. This results from tree
level annihilation of the quintuplets into $W$ or $Z$ bosons. The latter of
these arises from $\gamma Z$ and $ZZ$ final states, and is a simple
reweighting of the line cross section as

$\frac{\langle\sigma v\rangle_{ZZ}+\tfrac{1}{2}\langle\sigma
v\rangle_{Z\gamma}}{\langle\sigma
v\rangle_{\gamma\gamma}+\tfrac{1}{2}\langle\sigma v\rangle_{\gamma
Z}}=\frac{c_{\scriptscriptstyle W}^{2}}{s_{\scriptscriptstyle W}^{2}}.$ (85)

For the $W^{+}W^{-}$ final state, the tree level annihilation rate, with the
Sommerfeld effect included can be computed as,

$\displaystyle\langle\sigma v\rangle_{\scriptscriptstyle
WW}=\frac{\pi\alpha_{\scriptscriptstyle W}^{2}}{M_{\chi}^{2}}$
$\displaystyle\Big{[}18|s_{00}|^{2}+25|s_{0\pm}|^{2}+4|s_{0\pm\pm}|^{2}$ (86)
$\displaystyle+30\sqrt{2}{\rm Re}(s_{00}s_{0\pm}^{*})+12\sqrt{2}{\rm
Re}(s_{00}s_{0\pm\pm}^{*})+20{\rm Re}(s_{0\pm}s_{0\pm\pm}^{*})\Big{]}.$

As discussed for the case of the wino in Ref. Hryczuk:2011vi , higher order
corrections to this should not be appreciable, and so we do not include them.
We can then add the $W^{+}W^{-}$ final state to $dN/dE$ with weighting
$\langle\sigma v\rangle_{\scriptscriptstyle WW}/\langle\sigma v\rangle_{\rm
line}$ along with the $Z$ contribution. In each case, to determine the
spectrum of photons that result from these electroweak bosons we use PPPC4DMID
Cirelli:2010xx  with the electroweak corrections turned off, to avoid any
double counting of the endpoint contributions we computed. As seen in Fig. 8,
these contributions are important for $E_{\gamma}\ll M_{\chi}$.

The final contribution to the spectrum we consider is continuum photons
arising from the bound state decays. This contribution is not the main focus
of this work, but to get an estimate of its size and spectrum, we assume the
most common SM decay products are light quarks (equally weighted between
flavors), and employ the corresponding gamma-ray spectrum from PPPC4DMID. The
motivation for this choice is that the bound states will decay through their
couplings to (on-shell or off-shell) $W$ and $Z$ bosons, with the exact
channel depending on their $L$ and $S$ quantum numbers, and the gauge bosons
in turn decay dominantly to quarks due to their large associated degrees of
freedom. We weight the continuum spectrum by the ratio between the bound state
capture cross section and $\langle\sigma v\rangle_{\textrm{line}}$, similar to
how we weight the $Z$ and $W$ continuum components above. At the thermal mass,
this ratio is roughly 31, and the visible contribution can be seen in Fig. 8.
The contribution is dominated by the $Q=0$ $p\to s$ capture cross section
which sits near a Sommerfeld peak in this capture rate at 13.6 TeV (cf. Fig.
4). To highlight this, at the edges of the uncertainty band on the thermal
mass, 12.8 and 14.4 TeV, the equivalent ratio is reduced significantly, to
0.15 and 0.70, respectively. At 68.1 TeV the ratio is larger still – just over
471 – and is dominated by the $Q=0$ and $Q=1$ $p\to s$ capture rates.

Figure 9:  The full quintuplet spectrum for three different masses,
$M_{\chi}=14,\,16,$ and $18~{}{\rm TeV}$. As shown, the spectrum can change
significantly as a function of mass, a feature which does not arise for the
wino or higgsino.

Figure 8 highlights the various contributions to the spectrum, but does not
capture the variation of the spectrum as a function of mass. The variation can
be considerable, as shown in Fig. 9. From the definition, the line
contribution to this spectrum is fixed at $2\delta(E-M_{\chi})$ independent of
mass. What is not fixed, however, is the endpoint and continuum contributions,
which can vary significantly even for small changes in mass. (The bound state
contributions are not significant for the masses shown.) As shown in Ref.
Montanari:2022buj , such rapid variations can lead to sharp features in the
instrumental sensitivity to $\langle\sigma v\rangle_{\rm line}$, as the shape
of the DM signal being searched for varies rapidly with mass. These effects do
not occur for the wino or higgsino, where the spectra varies relatively
smoothly with mass (see Ref. Montanari:2022buj ).

The origin of this behavior seems to be interference between the different
Sommerfeld factors, associated with the distinct mass eigenstates for the
final annihilation: $\chi^{0}\chi^{0}$, $\chi^{+}\chi^{-}$, and
$\chi^{++}\chi^{--}$. These states have different branching ratios into the
various SM final states, and the positions of the resonance peaks can differ
between the interfering Sommerfeld factors. As we vary the mass, we can move
rapidly from the sharp turn off of one Sommerfeld peak to the sharp turn on of
another, and in doing so transition rapidly in the strength of the associated
endpoint and continuum contributions, as seen in Fig. 9.

One might ask why this behavior was not seen for the wino, which also has
multiple Sommerfeld factors that can interfere with each other. For the line
cross section, one might suspect that the issue is that only the
$\chi^{+}\chi^{-}$ state can annihilate to photons at tree-level; however, we
also do not see sharp quintuplet-like features in annihilation of wino DM to W
bosons, which is allowed at tree-level from both the $\chi^{0}\chi^{0}$ and
$\chi^{+}\chi^{-}$ states. More insight can be gained by working in the basis
of eigenstates of the potential at small $r$, rather than mass eigenstates;
this corresponds to the basis of potential eigenstates in the limit of
unbroken SU(2), as discussed in App. F. In the limit of unbroken SU(2), the
relevant potential for the quintuplet (coupling states with total charge zero
and even $L+S$) has two eigenstates that experience attractive interactions
and one eigenstate that experiences repulsive interactions. We expect that the
linear combination of Sommerfeld factors corresponding to the repulsed
eigenstate at small $r$ should be suppressed, as the SU(2) symmetry is
restored in the small-$r$ regime. We have confirmed numerically that this
suppression is quite pronounced, typically several orders of magnitude at the
velocities we consider. We would also expect a difference in the linear
combination of Sommerfeld factors corresponding to the two attracted
eigenstates, with the larger-magnitude eigenvalue yielding a larger Sommerfeld
enhancement, but this difference is much less dramatic, and so the two
attracted eigenstates still experience meaningful interference.

We thus attribute the sharp features observed in the quintuplet case to this
interference between the different attracted eigenstates (in the small-$r$
potential-dominated regime); its absence in the wino case is presumably
related to the fact that the wino has only one attracted eigenstate (e.g. Ref.
Asadi:2016ybp ). We thus expect this behavior (sharp variations in the
spectrum with mass) to be ubiquitous for larger SU(2) representations.

###  5.2 Uncertainty associated with the velocity distribution of dark matter

The complete initial-state wavefunctions naturally depend on the relative
velocity of the incoming DM particles, which in the discussion so far we have
simply set as $v=10^{-3}$. In this subsection we explore the systematic
uncertainties associated with our modeling of the velocity. We first discuss
the detailed dependence of the cross sections on the relative DM velocity, and
then explore the effects on our spectra of averaging over different plausible
velocity distributions. The effects of the long-range potential on the
wavefunction saturate when $v\lesssim m_{\scriptscriptstyle W}/M_{\chi}$,
which is true for halo velocities ($v\sim 10^{-3}$) for $M_{\chi}\lesssim
80~{}{\rm TeV}$; consequently, except near resonances or for very heavy DM, we
do not expect the Sommerfeld enhancement from the weak interactions to depend
sensitively on the velocity distribution. However, the bound-state formation
rate from $L>0$ partial-wave components of the initial state will have a non-
trivial velocity dependence even below this threshold. Furthermore, for the
quintuplet the thermal mass is only a factor of few below the saturation
threshold, and in systems with higher velocities than the Sun’s neighborhood –
such as galaxy clusters – both the direct annihilation cross section and the
bound-state formation rates are expected to depend sensitively on the
velocity.

Figure 10: The dependence of the cross sections on $v$ for four different
fixed masses. We show the case of direct endpoint annihilation (left, the
analogue of Fig. 11), and bound state capture (right, as in Fig. 12).

In Fig. 10 we show how annihilation and capture vary as a function of velocity
at four different masses. We observe that for a mass of $M_{\chi}=14$ TeV,
noticeable velocity dependence is present at $v\gtrsim 4\times 10^{-3}$. As we
discuss in depth in App. F, the oscillatory behavior observed at high
velocities can be understood in the limit of unbroken SU(2). This behavior
originates from interference between the different eigenvalues of the
potential. At low velocities, by contrast, SU(2) breaking effects are expected
to suppress the oscillations. For higher DM masses, where the velocity
dependence is relevant even for $v\lesssim 10^{-3}$, our previous annihilation
cross section plots should be taken as an illustrative estimate. A full
calculation at high mass would involve integrating the formulae given in this
paper over the true velocity distribution in the region of interest. The
oscillatory behavior of the cross section at high velocities means that
assuming a single velocity could in principle lead to large errors in this
case.

We now estimate the effect of averaging over the velocity distribution. The
characteristic scale of the DM velocity dispersion should be comparable to the
circular velocity of the visible matter, which in the vicinity of the Sun has
been measured to be $v_{\text{circ}}\simeq 240$ km/s 2016MNRAS.463.2623H .
Since the Milky Way’s rotation curve is roughly flat at the Sun’s location, we
expect the velocity dispersion to be of a similar order over much of the
Galaxy. However, close to the Galactic Center the DM velocity is not well-
known. In DM-only simulations the velocity dispersion falls as one approaches
the Galactic Center (e.g. Ref. Navarro:2008kc ) but simulations including
baryons have demonstrated the opposite behavior (e.g. Refs.
2010MNRAS.406..922T ; Board:2021bwj ; 2022MNRAS.513…55M ). Even at the Sun’s
location, the full DM velocity distribution is not well-understood: the
distribution is often treated as Maxwellian up to some escape velocity,
although this is only a crude approximation (e.g. Ref. Herzog-
Arbeitman:2017fte ). The escape velocity is determined to be $\sim$500 km/s at
the location of the solar system Necib:2021vxr ; kh2021 ; mfcEtAl2018
.171717Ref. mfcEtAl2018  finds that the escape velocity increases slightly
toward the Galactic Center. However, they only present results in to a radius
of around 5 kpc, where it is closer to 650 km/s. The precise value of this
cutoff is numerically unimportant for this work though, due to the exponential
suppression in the distribution (cf. Eq. (87)). In practice, we truncate the
particle velocity at 500 km/s, but the numerical difference between this and
2400 km/s is at most a part per million in the annihilation rate. Within the
Maxwellian approximation, the distribution is specified by that escape
velocity and the velocity dispersion, with the latter having a greater effect
on the annihilation rate.

For the Milky Way, we use the velocity dispersion values obtained in Ref.
Boddy:2018ike  for a variety of NFW-profiles. In particular, we take the
slowest and fastest velocities for locations interior to the solar system.
This gives a range $v_{\rm disp}\in[130,330]$ km/s. As a function of $v_{\rm
disp}$, the magnitude of the relative WIMP velocity is drawn from the
following 1D probability distribution,

$f(v)=\sqrt{\frac{27}{4\pi}}\,\frac{v^{2}}{v_{\rm disp}^{3}}e^{-3v^{2}/4v_{\rm
disp}^{2}}.$ (87)

Here $v_{\rm disp}$ is the RMS velocity for a single DM particle, which is
equal to the three-dimensional velocity dispersion $\sigma_{v,3d}$ defined in
Ref. Boddy:2018ike  by

$v_{\rm disp}^{2}=\sigma_{v,3d}^{2}=\frac{\int dv\,v^{4}f(r,v)}{\int
dv\,v^{2}f(r,v)}.$ (88)

Here $f(r,v)$ is the speed distribution for a single DM particle at a distance
$r$ from the Galactic Center.

Figure 11: An estimate of the range of uncertainty in our results associated
with the DM velocity distribution for the dominant direct annihilation. We use
two Maxwell distributions with $v_{\rm disp}$ – cf. Eq. (87) – at the extremal
values found by Ref. Boddy:2018ike . We divide their resulting $\langle\sigma
v\rangle$ by that of the simplified case of all quintuplets annihilating with
$v=10^{-3}$.

In Fig. 11, we plot the leading contribution to endpoint photon production,
direct annihilation from an $s$-wave initial state, for two different velocity
distributions, normalized to the simple assumption of all quintuplets having a
fixed $v=10^{-3}$. Except on resonances, the uncertainty is typically
negligible. We also see that off resonance, and particularly at lower masses,
the simple fixed-velocity assumption is a good approximation of either more
realistic model. The reason for this is simply that we are in the saturation
regime, as seen in Fig. 10. Therefore, we conclude that in general the fixed
velocity assumption is a good one at low masses, although at higher masses one
is generically underestimating the actual cross section, sometimes by more
than an order of magnitude. Accordingly, for an actual experimental analysis,
completeness would require an appropriate weighting of the cross section
according to the specific region of interest studied.

Figure 12: Capture rate by emission of a dipole photon from $p$-wave initial
state to $s$-wave quintuplet bound state for $Q=0$ states. The three choices
of velocity mirror those in Fig. 11. We see two different values of
$v_{\text{disp}}$ plotted, plus a $v_{\text{fixed}}=10^{-3}$ whose capture
rate we divide by.

For bound state capture, the off-resonance uncertainties are generally larger
than for direct annihilation, as anticipated. This is demonstrated in Fig. 12,
where we show $p$-to-$s$ capture, the leading single-photon dipole transition.
As we see by comparing the rates with those in Fig. 7 though, capture is
generally far subdominant to direct annihilation. In this channel, however,
the simple assumption of all DM having $v=10^{-3}$ generally provides a result
in the middle of the band given by the remaining two options. Even where it
does not, we see that it still provides a good approximation to the more
realistic velocity profiles. Accordingly, just as with direct annihilation, we
will take this value as a representative approximation of this subleading
contribution to endpoint photons when forecasting experimental sensitivity.

###  5.3 Estimating the experimental sensitivity to quintuplet DM

Finally we turn to an estimate of the experimental sensitivity to the
quintuplet DM hypothesis using the spectra we have computed. Using our
definition of the spectrum in Eq. (84), the average DM-generated flux an
instrument observes from a region of interest (ROI) of solid angle
$\Omega_{\rm ROI}$ is,

$\frac{d\Phi}{dE}=\frac{\langle\sigma v\rangle_{\rm line}}{8\pi
M_{\chi}^{2}}\left(\frac{dN}{dE}\right)_{\rm
smeared}\left(\frac{1}{\Omega_{\rm ROI}}\int
ds\,d\Omega\,\rho_{\chi}^{2}\right)\\!.$ (89)

As defined, the flux has units of $[{\rm counts}/{\rm cm}^{2}/{\rm s}/{\rm
TeV}/{\rm sr}]$ (for a detailed discussion of the units, see Ref.
Lisanti:2017qoz ). The final term in parentheses here is often referred to as
the $J$-factor, and is an integral over the DM density squared in the region
being observed. If the DM density $\rho_{\chi}$ was known exactly, then for a
model like the thermal quintuplet the flux is fully determined, as we have
computed both the cross section and spectrum (up to residual uncertainties
from higher order terms in the theory prediction, and velocity distribution as
discussed in the last subsection).

To test the quintuplet hypothesis, we need to compare the flux in Eq. (89) to
experimental measurements. For this, we will estimate the sensitivity of
H.E.S.S. to the endpoint photon signal using the “mock analysis” method
described in Ref. Baumgart:2017nsr . The approach in that work was to make use
of the publicly available H.E.S.S. data in Ref. HESS:2013rld , where a search
for $\chi\chi\to\gamma\gamma$ was performed using 112 hours of Galactic Center
observations taken by the instrument between $2004-2008$. This is a small
fraction of the observations H.E.S.S. has taken towards the galactic center.
As emphasized in Ref. Montanari:2022buj , the collaboration has already
collected roughly 800 hours of data in the region, and continues to collect
roughly 150 hours each year. In that sense, the dataset we consider represents
a small fraction of what is available. A further limitation of our approach is
that the analysis in Ref. HESS:2013rld  was a search purely for line photons,
and therefore adopted a flexible background model that absorbed all smooth
components. The analysis is therefore unsuitable to consider continuum
contributions, which can play an important role in these sorts of analyses, as
emphasized in, for instance, Ref. Rinchiuso:2020skh . Our rationale for
adopting this mock analysis, however, is that Ref. HESS:2013rld  provided
enough information that the full analysis they undertook can be performed
reconstructed, as was shown in Ref. Baumgart:2017nsr . Later on we will
provide a rough estimate of how sensitive more recent and upcoming analyses
could be.

Figure 13:  (Left) The estimated sensitivity of H.E.S.S. I using 112 hours of
galactic center (GC) data to the quintuplet. Assuming an Einasto profile, we
show sensitivity to $\langle\sigma v\rangle_{\rm line}$ as defined in Eq.
(84), which can then be compared to the equivalent theoretical prediction.
Across the entire mass range considered here, H.E.S.S. would be able to
exclude the quintuplet assuming an Einasto profile. (Right) If the DM profile
is cored as in Eq. (91), the core sizes that would be required to make a non-
observation consistent with our quintuplet predictions. At the thermal mass,
13.6 TeV, a core of 1 kpc is required, whereas at upper edge of the thermal
band, 14.4 TeV, the results would be consistent with a 0.5 kpc core.

For the mock analysis, we fit the data provided in Ref. HESS:2013rld  to a
combination of the flux in Eq. (89) and a parametric background model adopted
by the experiment—full details are provided in Ref. Baumgart:2017nsr . If we
have a prediction for $\rho_{\chi}$, this approach combined with our
prediction for the spectrum can be used to obtain an estimated sensitivity for
$\langle\sigma v\rangle_{\rm line}$. For this we adopt the Einasto profile
1965TrAlm…5…87E  used by the H.E.S.S. analysis (based on Ref. Pieri:2009je ),

$\rho_{\rm
Einasto}(r)\propto\exp\left[-\frac{2}{\alpha}\left(\left(\frac{r}{r_{s}}\right)^{\alpha}-1\right)\right]\\!,$
(90)

with $\alpha=0.17$, $r_{s}=20~{}{\rm kpc}$, and the normalization fixed to
$0.39~{}{\rm GeV/cm}^{3}$ at the solar radius, $r=8.5~{}{\rm kpc}$. The
resulting estimated constraint is shown in the left of Fig. 13.181818One can
find an early attempt to make projected $\gamma$-ray constraints from the
H.E.S.S. galactic center data in Ref. Cirelli:2015bda . Their Fig. 7 is
analogous to our Fig. 13. The earlier paper did not include LL (or NLL)
resummation, nor the NLO corrections to the electroweak potential. For similar
Einasto parameters, their bounds on the line cross section are about an order
of magnitude weaker than ours (however, note that our predicted line cross
section is also smaller). We emphasize, that even though this is a constraint
on $\langle\sigma v\rangle$ it is not based solely on the line prediction—the
results are based on the estimated detectability of the entire spectrum
resulting from line, endpoint, and bound state photons. For the mass range
considered, the bound state contribution is negligible (cf. Fig. 7), however
the endpoint is not: at 13.6 TeV, it enhances the sensitivity by a factor of
1.9. As mentioned above, by default, the results do not include either
continuum contribution considered in Fig. 8. The motivation for this is the
particular background model adopted in Ref. HESS:2013rld  was designed solely
to search for a narrow line feature, and there can be considerable degeneracy
with the continuum emission (see the discussion in Ref. Baumgart:2017nsr ).
Nevertheless, we have tested adding the continuum emission from direct
annihilation to $W$ and $Z$ final states, and found for most masses it has no
impact on the estimated limits, although there is a slight fluctuation around
the thermal mass, which increases the sensitivity by $\simeq$$20\%$. We also
note that there are several locations, such as just below $M_{\chi}=3~{}{\rm
TeV}$ and just above the thermal mass where there is a larger theoretical
error. This results from the sharp variations in the spectra observed in Fig.
9. In fact, the sensitivity to these features is reduced by the insensitivity
of the background model used in this work to smooth features. These results
can also be compared to Ref. Montanari:2022buj , where an alternate H.E.S.S.
analysis is performed using the spectra from this work, and much sharper
variations are observed in their sensitivity. (We note that the results of
that work made use of the LO Sommerfeld potential calculations, not the NLO
results we have adopted here.)

The results of the mock analysis suggest that for the central value of the
thermal mass, even 112 hours of H.E.S.S. data can exclude the thermal
prediction by a factor of 10. Nevertheless, this varies considerably across
the uncertainty band on the thermal mass: at 12.8 TeV, the exclusion factor is
55, whereas it is only 4 at 14.4 TeV. The sharp variation is a result of the
thermal window sitting near a Sommerfeld resonance, as shown in Fig. 13. We
emphasize once more that although we use the NLO potential in our
computations, the thermal mass range was computed with the LO potential, and
given the sensitivity of our findings to the exact mass, computing the thermal
mass at NLO will be important for narrowing the fate of the thermal
quintuplet. If we relax the thermal cosmology assumption and consider a
broader range of masses, we see the quintuplet is excluded across the full
$0.5-20~{}{\rm TeV}$ mass range. Yet this statement is contingent on the form
of $\rho_{\chi}$ adopted, about which there is considerable uncertainty. In
particular, the density profile may flatten toward the inner Galaxy. As the
annihilation signal is sensitive to $\rho_{\chi}^{2}$ flattening the profile
has a marked impact on the flux, making this one of the dominant uncertainties
in $\rho_{\chi}$ for our purposes. We parameterize a possible flattening of
the profile by replacing the Einasto density profile with a constant value for
Galactocentric distances $r<r_{c}$, where we will refer to $r_{c}$ as the
“core size”:191919To be explicit, we fix the normalization of the Einasto
profile in Eq. (90) to the solar radius before we impose the core restriction
of Eq. (91). This implies that for a core size larger than the solar radius,
the profile will predict less than 0.39 GeV/cm3 at our location. Of course,
the more significant concern is that such a large core is not consistent with
observations, as discussed in the text, and should solely be viewed as a proxy
for how much the $J$-factor needs to be reduced.

$\rho(r)=\begin{cases}\rho_{\rm Einasto}(r)&r>r_{c},\\\ \rho_{\rm
Einasto}(r_{c})&r<r_{c}.\end{cases}$ (91)

We can then ask what choice of $r_{c}$ would raise the estimated constraint on
$\langle\sigma v\rangle_{\text{line}}$ above the theoretical prediction. This
is plotted in the right-hand panel of Fig. 13, both employing our full
endpoint spectrum and in the case where we (incorrectly) use only the line
cross section in setting the bounds. To provide an estimate for what core
sizes are consistent with data, we note that simulations of Milky Way like
galaxies can generate ${\cal O}(1~{}{\rm kpc})$ cores Chan:2015tna , however,
measurements of stars in the Bulge seem to disfavor $r_{c}\gtrsim 2~{}{\rm
kpc}$ 2015MNRAS.448..713P ; Hooper:2016ggc . At the lower end of the thermal
mass range, 12.8 TeV, the thermal Quintuplet would already be in tension with
this, requiring a 2.8 kpc core. At the central (13.6 TeV) and upper (14.4 TeV)
end, however, the required core size is 1.0 and 0.5 kpc, respectively, and
therefore not obviously in tension with our mock analysis. A more recent study
claims evidence for a few kpc core that could potentially saturate the earlier
limits Ou:2023adg . If confirmed, this could suppress the $J$-factor by nearly
an order of magnitude. That would challenge the indirect-detection community
to set aggressive limits, but as we estimate, even cores of this size are in
reach of CTA and possibly H.E.S.S (cf. Fig. 14).

As already mentioned, the mock analysis we consider makes use of only a very
small amount of existing data, and with that we forecast a sensitivity to
$\langle\sigma v\rangle_{\rm line}\simeq 8.5\times 10^{-27}~{}{\rm cm^{3}/s}$
at the central thermal mass. (The error on this value due to uncertainty on
the spectrum from our NLL calculation is less than 1%, far smaller than the
variation across the thermal mass range, which is closer to 10%.) Using 500
hours of H.E.S.S. data and an identical Einasto profile, Ref.
Montanari:2022buj  forecast a sensitivity at the thermal mass of
$\simeq$$9.3\times 10^{-28}~{}{\rm cm^{3}/s}$, almost a factor of ten better
than used here. This is significantly more than the naive $\sqrt{5}$ the
additional data would suggest, which can be primarily attributed to the fact
that work used H.E.S.S. II observations, whereas we adopt the sensitivity from
H.E.S.S. I, combined with a different analysis used in that work. With such a
sensitivity, we would require a core size slightly larger than $3.5~{}{\rm
kpc}$ to save the thermal quintuplet, which would be in tension with
observations. Nevertheless, repeating the process at 14.4 TeV, the required
core size would only be $1.6~{}{\rm kpc}$, and therefore not yet clearly
excluded. We can also give a crude estimate for the sensitivity the upcoming
CTA could have for the quintuplet. Although no dedicated forecast for the
quintuplet has been performed using the spectra in our work, we can estimate
the improved sensitivity as follows. Reference Baumgart:2018yed  performed an
identical mock analysis to ours for the NLL wino spectrum, estimating
sensitivity at $M_{\chi}=13.6~{}{\rm TeV}$ of $\langle\sigma v\rangle_{\rm
line}\simeq 8\times 10^{-27}~{}{\rm cm^{3}/s}$, slightly stronger than the
sensitivity to the quintuplet. Using the identical NLL spectrum, Ref.
Rinchiuso:2020skh  then estimated that with 500 hours of data, CTA could reach
$\simeq 1\times 10^{-28}~{}{\rm cm^{3}/s}$, a factor of eighty improvement.
Assuming the same improvement for the quintuplet, CTA would be sensitive to
$\langle\sigma v\rangle_{\rm line}\simeq 1.1\times 10^{-28}~{}{\rm cm^{3}/s}$,
excluding the thermal value by a factor of eight hundred. To not have seen the
thermal quintuplet, we would need to core $\rho_{\chi}$ out to almost
$8.6~{}{\rm kpc}$ – beyond the solar radius – which is simply inconsistent
with observations. Even at the upper end of the mass range, a core of
$6.4~{}{\rm kpc}$ would be required. In this sense, CTA would provide the
definitive word on the whether the thermal quintuplet is the DM of our
Universe.

Figure 14: An estimate for the required core size of the Einasto profile,
given no preference for a signal is seen in our H.E.S.S. mock analysis, and if
no signal emerges in an analysis of the full H.E.S.S. dataset, or at CTA. The
central values correspond to 13.6 TeV, the central thermal mass, whereas the
upper and lower error bars correspond to 12.8 and 14.4 TeV, the edges of the
thermal mass window. The dashed corresponds to the rough core size constraint
from Ref. Hooper:2016ggc . Our results suggest that H.E.S.S. can already
considerably test the quintuplet, with the final word likely being left to
CTA.

These results are summarized in Fig. 14, where the point shows the core size
required for the central thermal mass, 13.6 TeV, whereas the upper and lower
error bars corresponds to the lower and upper ends of the thermal mass range.
We also show the core size disfavored by the analysis in Ref. Hooper:2016ggc .
The figure summarizes the conclusion reached above: CTA will seemingly have
the final word of the thermal quintuplet, however if a full analysis of
H.E.S.S. data sees no sign of the signal, the model would already begin to be
disfavored. There are two important caveats to this conclusion. The first is
that our findings are based off an extrapolation from a mock analysis of
H.E.S.S. I data, and are no substitution for a full analysis or projection
using the present and forecast H.E.S.S. and CTA instrumental responses. To
give one example of what could change, an analysis that accounts for the
continuum emission could be able to even more strongly test the quintuplet.
Secondly, the range of masses we have considered is the thermal mass window of
$13.6\pm 0.8~{}{\rm TeV}$ that was determined using the LO potential, whereas
the remainder of our calculations use the NLO results, as emphasized several
times already. Updating the thermal mass using the NLO potentials will be
important. To give a sense for the impact this could have, repeating the
analysis in Fig. 14 for 13.6 TeV using the LO potential in our calculations,
the core size would change from 1.0, 3.5, and 8.5 kpc, to 0.7, 2.5, and 7.1
kpc.

##  6 Conclusions

For all the vastness of the DM parameter space, a thermal WIMP has remained a
constant focus for decades. Minimal DM is an exemplar of thermal DM, and
through indirect detection, many of the associated models are on the verge of
being detected or firmly excluded, as we have shown for the quintuplet in the
present work. Either way, these are important times in the search for DM.

With this in mind, the present work has computed the quintuplet annihilation
spectrum to NLL accuracy, and established the formalism to straightforwardly
extend this to higher odd SU(2) representations. We plot the spectrum along
with projected limits from a simple extension of the H.E.S.S. I analysis in
Fig. 13. In doing so, we have demonstrated the power of the EFT of Heavy DM,
and also extended this formalism to include the contribution from the rich set
of bound states the model contains. While the bound states can make a
significant contribution to the continuum photon emission, their impact on the
number of photons with $E_{\gamma}\sim M_{\chi}$ is minimal, except at
isolated masses. As seen in earlier studies of the wino, the same cannot be
said for endpoint photons from direct annihilation, which again provide an
${\cal O}(1)$ correction to the line signal seen in IACTs.

Taken together, we estimated that with our spectra, H.E.S.S. should almost be
able to probe the entire allowed range for the quintuplet once uncertainties
on the DM density in the inner galaxy are accounted for. Performing this
analysis using the existing data, and the soon-to-be-collected data with CTA
will be critical (cf. Fig. 14). The use of background models which enhance
sensitivity to smooth features such as the continuum as well as the full
contribution from the endpoint-photon spectrum can provide an additional piece
of experimental leverage beyond conventional line searches.

On the theory side, the thermal abundance should be recomputed using NLO
potentials, as the sensitivity to the quintuplet depends strongly on what end
of the predicted thermal mass range one sits. Finally, it will be interesting
to extend the techniques in this work to additional representations, such as a
$\mathbf{7}$ of SU(2), where we expect key features of the quintuplet such as
the strong variation in the spectrum as a function of mass, to appear and be
even more pronounced.

###### Acknowledgements.

The work we have presented benefited from useful discussions with Tobias
Binder, Marco Cirelli, Tongyan Lin, Alessandro Montanari, Emmanuel Moulin,
Nikhil Raghuram, and Diego Redigolo. MB is supported by the DOE (HEP) Award
DE-SC0019470. VV is supported by startup funds from the University of South
Dakota. TRS’ work is supported by the Simons Foundation (Grant Number 929255,
T.R.S), by the National Science Foundation under Cooperative Agreement
PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental
Interactions), and by the U.S. Department of Energy, Office of Science, Office
of High Energy Physics of U.S. Department of Energy under grant Contract
Number DE-SC0012567. TRS thanks the Aspen Center for Physics, which is
supported by National Science Foundation grant PHY-2210452, for hospitality
during the completion of this work.

## Appendix

##  Appendix A Quintuplet Dark Matter: A Brief Review

Here we provide a brief review of quintuplet DM, also referred to as 5-plet
electroweak DM. In particular, we firstly outline the relevant group theory
necessary to specify the interactions used for the calculations in the main
text. We will then review phenomenological aspects of the model beyond
indirect detection.

###  A.1 Interactions

Quintuplet DM consists of adding to the SM five Majorana fermions that
transform together in the $\mathbf{5}$ representation of SU(2), and as a
singlet under the remaining SM forces. Above the electroweak symmetry breaking
scale, we collect the five fields into a multiplet
$\chi=(\chi^{1},\,\ldots,\,\chi^{5})^{T}$, in terms of which the DM Lagrangian
takes the following form

$\displaystyle\mathcal{L}_{\scriptscriptstyle\textrm{DM}}=\frac{1}{2}\bar{\chi}(i\not{D}-M_{\chi})\chi=\frac{1}{2}\bar{\chi}([i\not{\partial}-M_{\chi}]\mathbb{1}+g_{\scriptscriptstyle
W}T_{\mathbf{5}}^{a}\not{W}^{a})\chi.$ (92)

In the final expression, the first two contributions represent the kinetic
terms for each of the fields, which are diagonal amongst the multiplets as
indicated by $\mathbb{1}$. The final term describes the interaction between
the additional fields and the SM electroweak bosons. Importantly, we emphasize
that the interaction strength is specified by the SM SU(2) gauge coupling,
$g_{\scriptscriptstyle W}$, and is not a free parameter. Instead, $M_{\chi}$
is the unique free parameter in the theory. If we assume a conventional
thermal origin for the quintuplet, then even the mass can be fixed by the
observed relic density to $M_{\chi}=13.6\pm 0.8~{}{\rm TeV}$ Mitridate:2017izz
; Bottaro:2021snn .202020As emphasized above, this value was computed with the
LO electroweak potential. Redoing the analysis with the NLO potential would
reduce an important theoretical uncertainty. The 13.6 TeV quintuplet is
therefore a zero parameter DM model.

After electroweak symmetry breaking, the five Majorana fermions rearrange
themselves into three mass eigenstates: a Majorana neutral fermion $\chi^{0}$,
and two charged Dirac fermions $\chi^{+}$ and $\chi^{++}$. At leading order,
each of these states maintains a mass of $M_{\chi}$. However, radiative
corrections to the charged states break this degeneracy, raising the masses of
the charged fermions, singling out $\chi^{0}$ as the lightest state and the DM
candidate. These corrections have been computed, and in detail
$\delta_{0}=M_{\chi^{+}}-M_{\chi^{0}}\simeq 164~{}{\rm MeV}$, and
$\delta_{+}=M_{\chi^{++}}-M_{\chi^{+}}=3\delta_{0}$ Cirelli:2005uq ;
Ibe:2012sx . For most aspects of our calculations these mass splittings will
be irrelevant and we will take $\delta_{0}\simeq\delta_{+}\simeq 0$. However,
we do include them in the electroweak potential used to compute Sommerfeld
factors, and scattering- & bound-state wavefunctions. This is done by adding
$2\delta_{0}$ to the diagonal term in the potential matrix correspondingto the
$\chi^{+}\chi^{-}$ component of any state, $8\delta_{0}$ to the diagonal term
corresponding to the $\chi^{++}\chi^{--}$ component, and similar, appropriate
shifts for diagonal elements in the potential matrix corresponding to
components of $Q\neq 0$ states (the shift is given by the difference between
the rest mass of the state constituents and $\chi^{0}\chi^{0}$).

We now turn to the interaction term in Eq. (92),
$\tfrac{1}{2}g_{\scriptscriptstyle
W}\bar{\chi}T_{\mathbf{5}}^{a}\not{W}^{a}\chi$. Here $a=1,2,3$ indexes the
electroweak gauge bosons, which transform together in an adjoint of SU(2). In
the broken theory, these are mapped to the charge and mass eigenstates in the
usual way,

$\displaystyle
W^{1}_{\mu}=\frac{1}{\sqrt{2}}(W^{+}_{\mu}+W^{-}_{\mu}),\;\;W^{2}_{\mu}=\frac{i}{\sqrt{2}}(W^{+}_{\mu}-W^{-}_{\mu}),\;\;W^{3}_{\mu}=s_{\scriptscriptstyle
W}A_{\mu}+c_{\scriptscriptstyle W}Z_{\mu}.$ (93)

The only part of Eq. (92) that remains undetermined is $T_{\mathbf{5}}^{a}$,
the three generators of SU(2) in the quintuplet representation. A convenient
basis in which to specify $T_{\mathbf{5}}^{a}$ is the basis of charged states
discussed above, where the DM can be cleanly identified. We can determine the
charged states through their couplings to the bosons in Eq. (93). In
particular, as $A_{\mu}$ couples to charge, we can read off the charges of the
states as soon as we know $T^{3}$. It will also be convenient to introduce

$T^{\pm}=\frac{1}{\sqrt{2}}(T^{1}\pm iT^{2}),$ (94)

in terms of which $T^{a}W^{a}=T^{+}W^{+}+W^{-}T^{-}+W^{3}T^{3}$, independent
of representation.

Before evaluating the generators for the quintuplet, let us review how the
argument proceeds for the simpler case of the wino—a triplet or $\mathbf{3}$
of SU(2). In this case, it is conventional to exploit the fact that the
generators are given by the structure constants of SU(2), so that

$\bar{\chi}T_{\mathbf{3}}^{a}\gamma^{\mu}\chi=\bar{\chi}_{b}(T_{\mathbf{3}}^{a})_{bc}\gamma^{\mu}\chi_{c}=-i\epsilon_{abc}\bar{\chi}_{b}\gamma^{\mu}\chi_{c}.$
(95)

This approach depended on the fact we already knew a representation of the
generators for the adjoint, and so does not simply generalise to larger
representations. However, we can derive a representation more systematically
as follows. Recall that one representation of $\mathbf{n}$ in SU(2) is an
$n-1$ index symmetric tensor, with each index transforming in the fundamental.
In this representation we denote the adjoint as $\chi^{ij}$, and the
quintuplet as $\chi^{ijkl}$, where the indices take values 1 and 2. Beginning
with the wino, $\chi^{ij}$ has three unique components, which we can embed
into a vector as212121The $\sqrt{2}$ ensures that $\chi$ as it appears in
$\bar{\chi}\not{\partial}\chi/2=\bar{\chi}^{ij}\not{\partial}\chi^{ij}/2$ is
canonically normalized. An identical argument explains the coefficients in Eq.
(100).

$\chi=\begin{pmatrix}\chi^{1}\\\ \chi^{2}\\\
\chi^{3}\end{pmatrix}=\begin{pmatrix}\chi^{11}\\\ \sqrt{2}\chi^{12}\\\
\chi^{22}\end{pmatrix}\\!.$ (96)

We can use this representation to explicitly construct $T^{a}_{\mathbf{3}}$ as
follows. The key is that we know exactly how the generators act on $\chi^{ij}$
as each index transforms in the fundamental. Accordingly,

$T^{a}(\chi^{ij})=[(T^{a}_{F})^{i}_{k}\delta^{j}_{l}+\delta^{i}_{k}(T^{a}_{F})^{j}_{n}]\chi^{kl},$
(97)

where $T^{a}_{F}=\sigma^{a}/2$, with $\sigma^{a}$ the Pauli matrices. Consider
a generic infinitesimal transformation, $U=1+iu$, with $u=u_{a}T^{a}$. If we
take $u^{a}=(0,0,\kappa)$ with $\kappa\ll 1$, then the components of $\chi$
transform as

$\displaystyle\delta\chi^{1}=i\kappa\chi^{1},\;\;\delta\chi^{2}=0,\;\;\delta\chi^{3}=i\kappa\chi^{3}\\!.$
(98)

From this, we can read off the action of an infinitesimal $\mathbf{U}$ on
$\chi$, and hence infer that we must have

$T^{3}_{\mathbf{3}}={\rm diag}(+1,\,0,\,-1).$ (99)

We can now identify $\chi^{1}$, $\chi^{2}$, and $\chi^{3}$ as having charges
$+1$, $0$, and $-1$, yielding the expected spectrum in the broken
phase.222222The representation in the charge basis is not unique. For
instance, the transformation $\chi^{1,2}\to e^{\pm i\phi}\chi^{1,2}$ leaves
the charge assignments unchanged, but will introduce a phase into the off-
diagonal $W^{\pm}$ couplings. The same will be true for the off-diagonal
quintuplet couplings. The remaining components of $T^{a}_{\mathbf{3}}$ can be
derived identically.

This approach readily generalizes to the quintuplet. The five unique
components of $\chi^{ijkl}$ can be embedded into a vector as follows,

$\chi=\begin{pmatrix}\chi^{1}\\\ \chi^{2}\\\ \chi^{3}\\\ \chi^{4}\\\
\chi^{5}\end{pmatrix}=\begin{pmatrix}\chi^{1111}\\\ 2\chi^{1112}\\\
\sqrt{6}\chi^{1122}\\\ 2\chi^{1222}\\\
\chi^{2222}\end{pmatrix}=\begin{pmatrix}\chi^{++}\\\ \chi^{+}\\\ \chi^{0}\\\
\chi^{-}\\\ \chi^{--}\end{pmatrix}\\!.$ (100)

To justify the charge assignments, we repeat the above analysis to find

$T^{3}_{\mathbf{5}}={\rm diag}(+2,\,+1,\,0,\,-1,\,-2).$ (101)

Further, we can compute

$\displaystyle T^{+}_{\mathbf{5}}=\begin{pmatrix}0&\sqrt{2}&0&0&0\\\
0&0&\sqrt{3}&0&0\\\ 0&0&0&\sqrt{3}&0\\\ 0&0&0&0&\sqrt{2}\\\
0&0&0&0&0\end{pmatrix}\\!,$ (102)

and $T^{-}_{\mathbf{5}}=(T^{+}_{\mathbf{5}})^{T}$. In the main body we will
exclusively work in the charged basis of Eq. (100), using the form of
$T_{\mathbf{5}}^{a}$ as above whenever necessary.

###  A.2 Phenomenology

We end this section with a brief description of quintuplet phenomenology
beyond indirect detection. As already noted, the case where $M_{\chi}\simeq
13.6$ TeV is particularly appealing, as this is the mass singled out from a
conventional cosmology via a WIMP-miracle-style argument. The interactions of
the quintuplets, reviewed in the previous subsection, are sufficient to keep
it in thermal equilibrium in the early Universe. As the Universe cools,
eventually the quintuplet undergoes a conventional freeze-out. When this
occurs dictates the final abundance, and matching to the observed DM density
fixes the single parameter of the theory, $M_{\chi}$. For larger (smaller)
values of $M_{\chi}$ than 13.6 TeV, one naively over (under) produces the
observed DM density. Nevertheless, it is entirely possible that there are
effects in the early Universe that modify this simple picture. The presence of
additional beyond-the-SM states that either decay to the SM or directly to the
quintuplet can dilute or increase its abundance, respectively, making a wider
mass range viable. Further, for $M_{\chi}\leq 13.6$ TeV, even with no
additional states, the quintuplet would represent a well motivated (and
predictable) fraction of DM. For these reasons, in the main body we considered
a wide range of quintuplet masses, although we emphasize once more that the
scenario where $M_{\chi}\simeq 13.6$ TeV is compelling.

As with any DM candidate, one can also consider searching for the quintuplet
with either direct detection or at a collider. For direct detection, as the
quintuplet carries no U(1) hypercharge, it does not couple to the $Z$ boson at
tree level. Nevertheless, couplings to SM nucleons can arise at loop level.
The spin-independent cross-section is $\sigma\simeq(1.0\pm 0.3)\times
10^{-46}~{}{\rm cm}^{2}$ Bottaro:2021snn  (see also Ref. Hisano:2015rsa ;
Chen:2023bwg ), beyond the reach of current searches at 13.6 TeV, and even
next generation instruments such as LZ Mount:2017qzi . Nevertheless, the
cross-section is a factor of $\sim 4$ above the neutrino floor, and is in
reach of generation-3 instruments such as DARWIN DARWIN:2016hyl .

Detection at colliders is similarly challenging, but also potentially within
reach of future instruments. Existing LHC searches reach masses of $\sim
270~{}\text{GeV}$; even the high-luminosity dataset will only reach to $\sim
520~{}\text{GeV}$ Ostdiek:2015aga . Even future hadron colliders are unlikely
to reach the thermal mass. A future 100 TeV hadron collider will just be able
to reach the thermal masses for canonical neutralinos such as the higgsino and
wino Cirelli:2014dsa ; Low:2014cba . Given these two candidates have
significantly lower thermal masses of 1 and 2.9 TeV respectively Hisano:2006nn
; Cirelli:2007xd ; Hryczuk:2010zi ; Beneke:2016ync , the prospects of probing
the 13.6 TeV quintuplet appear discouraging. Future muon colliders operating
at lower center-of-mass energies could reach the quintuplet, although they
would still need to obtain $\sqrt{s}\simeq 35~{}{\rm TeV}$ Bottaro:2021snn
(see also Refs. Han:2020uak ; Han:2022ubw ). Taken together, in the short term
indirect detection remains the most likely avenue for probing quintuplet DM.

##  Appendix B Operators for higher-$L$ Bound State annihilation

As argued in Sec. 4.1, the higher-$L$ bound states preferentially decay to the
deeper bound states with lower $L$ instead of directly annihilating to SM
particles. Nevertheless, for completeness here we provide the complete set of
relevant operators up to ${\cal O}(v)$.

We consider the structure of the operators which support up to $p$-wave bound
states. Thus we need to keep operators (at the amplitude level) suppressed by
at most one power of the DM 3-momentum; the $\mathcal{O}(v^{0})$ operators
will support both the direct annihilation as well as $s$-wave bound state
annihilation. By matching to the full tree level amplitude, we can obtain the
structure that supports $S=0,\,L=1$ states,232323Throughout this appendix we
will work with 4-component DM fields and not reduce it to 2 components as was
done for the $L=S=0$ operator in the main text.

$\mathcal{O}_{1}=\mathbf{v}_{\chi}\cdot\mathbf{n}\left(\bar{\chi}\Big{[}T^{a},T^{b}\Big{]}\gamma^{0}\gamma^{5}\chi\right)i\epsilon^{ijk}(n-\bar{n})^{k}{\cal
B}_{\perp n}^{i,a}{\cal B}_{\perp\bar{n}}^{j,b},$ (103)

where it is understood that the subscript $\chi$ on $\mathbf{v}_{\chi}$
indicates that the velocity vector is the velocity of the state
created/annihilated by the $\chi$ field (as opposed to the $\bar{\chi}$
field). As is evident, this operator supports a bound state with $L=1$ and
$S=0$. Likewise we can write down the next set of operators,

$\displaystyle\mathcal{O}_{2}$
$\displaystyle=\mathbf{v}_{\chi}\cdot\mathbf{n}\left(\bar{\chi}\Big{\\{}T^{a},T^{b}\Big{\\}}\gamma^{i}\chi\right)(n-\bar{n})^{i}{\cal
B}_{\perp n}^{\mu,a}{\cal B}^{b}_{\perp\bar{n}\mu},$ (104)
$\displaystyle\mathcal{O}_{3}$ $\displaystyle=\mathbf{v}_{\chi}\cdot{\cal
B}^{a}_{\perp n}{\cal
B}^{b}_{\perp\bar{n}\mu}\left(\bar{\chi}T^{b}T^{a}\gamma^{\mu\perp}\chi\right)\\!,$
$\displaystyle\mathcal{O}_{4}$ $\displaystyle=\mathbf{v}_{\chi}\cdot{\cal
B}^{b}_{\perp\bar{n}}{\cal B}^{a}_{\perp
n\mu}\left(\bar{\chi}T^{a}T^{b}\gamma^{\mu\perp}\chi\right)\\!.$

From their Dirac structure, it is clear that these operators support $L=1,S=1$
bound states. There is also another operator which supports an $L+S$ odd bound
state, specifically the $L=1$, $S=0$ bound state which arises out of a
correction to an ultra-soft gauge boson emission off the heavy DM particles.
The details of this operator are involved and hence are separately given
further below.242424A complete anaysis and resummation of this channel at NLL
would require us to do a full two loop computation in order to recover the
anomalous dimension as well matching to stage II of the EFT. We leave this
work for the future.

We note that decays of SU(2)-singlet bound states with $L+S$ odd into two
(transverse) gauge bosons are constrained by charge conjugation invariance and
parity. This is because the underlying Lagrangian in Eq. (92) which couples
the electroweak gauge sector and a Majorana quintuplet fermion is $C$ and $P$
invariant. A fermion-antifermion bound state has $C$ eigenvalue $(-1)^{L+S}$
and $P$ eigenvalue $(-1)^{L+1}$. A $\gamma\gamma$, $\gamma Z$, or $ZZ$ final
state has $C$ eigenvalue $+1$, thereby forbidding an $L+S$ odd bound state
decay into them by $C$ alone. Decays to the $W^{+}W^{-}$ final state are
allowed regardles of the $L,\,S$ quantum numbers of the bound state.252525One
might think the Landau-Yang theorem also forbids decay from $L+S$-odd initial
states into two bosons at all orders, but the application of the theorem to
non-Abelian theories is subtle; there is a generalized Landau-Yang theorem
that holds so long as the decay products are in a color-singlet state
Beenakker:2015mra , but this is not the case for the decay of $L+S$-odd
states. In Ref. Asadi:2016ybp , one can see explicitly that for wino-onium,
decays to $W^{+}W^{-}$ are allowed for all combinations of initial-state $L$
and $S$. In all cases, decays to longitudinal gauge bosons are potentially
allowed, as well.

Finally, let us return to the additional operator which is sub-leading in
velocity that we can write down by looking at the emission of another gauge
boson which usually contributes to the $Y_{v}$ Wilson line. We do not get any
such correction from the $Y_{n}$ or the $Y_{\bar{n}}$ Wilson line since we do
not wish to consider sub-leading terms in the SCET power counting parameter.

We start again with our ${\cal O}(v^{0})$ operator at the amplitude level
before we dress it with soft Wilson lines

$\displaystyle\bar{\chi}\\{T_{\chi}^{a},T_{\chi}^{b}\\}\Gamma\chi{\cal
B}_{\perp n}^{a}{\cal B}^{b}_{\perp\bar{n}},$ (105)

where $\Gamma$ is an arbitrary Lorentz structure. Consider the emission of an
SCET ultra-soft gauge boson of momentum k off the initial $\chi$ or
$\bar{\chi}$ particle. First looking at the $\chi$ particle,

$\displaystyle
u_{e}(\tilde{p})+ig\frac{i(\not{p}-\not{k}+M_{\chi})}{(p-k)^{2}-M_{\chi}^{2}+i\epsilon}\gamma^{\mu}\epsilon_{\mu}^{a}(T_{\chi}^{a})_{ce}u_{e}(p)$
(106)

$\tilde{p}$, p is the momentum of the $\chi$ particle. Let us now expand this
result to $\mathcal{O}(v)$ (except inside the spinor). We see that apart from
the usual Wilson line contribution, we also get two other relevant terms,

$\displaystyle u_{e}(\tilde{p})-g_{\scriptscriptstyle
W}(T_{\chi}^{a})_{ce}\left(-\frac{\epsilon_{0}^{a}}{k_{0}}+\frac{\mathbf{v}\cdot\mathbf{\epsilon}^{a}(k)}{k_{0}}-\frac{\epsilon_{0}^{a}(k)\mathbf{v}\cdot\mathbf{k}}{k_{0}^{2}}\right)u_{e}(M_{\chi}).$
(107)

If we sum the infinite series of gauge bosons with one of the propagators
expanded out to $\mathcal{O}(v)$, we then have a structure

$\displaystyle\chi\rightarrow Y_{v}\mathbf{v}\cdot\mathbf{B_{s}}\chi$ (108)

where we have defined the following Hermitian operator,

$\displaystyle\mathbf{B}_{s}=\frac{Y^{\dagger}_{v}\left(\mathbf{\mathcal{P}}-g_{\scriptscriptstyle
W}\mathbf{A}_{s}\right)Y_{v}}{v\cdot\mathcal{P}},$ (109)

where $\mathcal{P}$ is the momentum label operator. If this term is included
in a larger expression, the $v\cdot\mathcal{P}$ factor in the denominator only
acts on the terms in the numerator while the one in the numerators acts only
on the $Y_{v}$ Wilson line to the right. We can now combine this with the
emissions off $\bar{\chi}$ to give us an effective operator, now dressed with
soft Wilson lines

$\displaystyle\bar{\chi}\Big{\\{}Y_{v}^{\dagger}\\{T^{a},T^{b}\\}\Gamma
Y_{v},\mathbf{v}\cdot\mathbf{B}_{s}\Big{\\}}\chi{\cal B}_{\perp
n}^{a^{\prime}}{\cal
B}^{b^{\prime}}_{\perp\bar{n}}Y_{n}^{aa^{\prime}}Y_{\bar{n}}^{bb^{\prime}}.$
(110)

We can once again use our Wilson line identity to write this in terms of our
usual soft function

$\displaystyle\bar{\chi}\Big{\\{}\\{T^{a},T^{b}\\}\Gamma,\mathbf{v}\cdot\mathbf{B}_{s}\Big{\\}}\chi{\cal
B}_{\perp n}^{c}{\cal B}^{d}_{\perp\bar{n}}Y^{abcd}.$ (111)

We can then formally define a new object $\mathbfcal{B}_{s}$ to explicitly
separate out all the soft fields from the heavy DM fields

$\displaystyle\mathbfcal{B}_{s}^{a}T^{a}=\mathbf{B}_{s}.$ (112)

Then we see that

$\displaystyle\mathbfcal{B}_{s}^{a}=\mathrm{Tr}[\mathbf{B}_{s}T^{a}],$ (113)

where we have used $\mathrm{Tr}[T^{a}T^{b}]=\delta^{ab}$. Our operator now
becomes

$\displaystyle\bar{\chi}\Big{\\{}\\{T^{a},T^{b}\\}\Gamma,T^{e}\Big{\\}}\chi\,{\cal
B}_{\perp n}^{c}{\cal
B}^{d}_{\perp\bar{n}}Y^{abcd}\mathbf{v}\cdot\mathbfcal{B}_{s}^{e}.$ (114)

In summary, we have an additional soft function which is explicitly suppressed
in velocity. The Dirac structure for this operator is $\gamma^{0}\gamma^{5}$
so that this is the only operator that supports an L+S odd bound state.
However, it is clear that the soft operator only begins at one loop and hence
we are always forced into at least a 3 gauge boson final state.

At the amplitude squared level, we simply have a square of this operator and
there will no interference terms with other bound state operators since all
other operators support an L+S even bound state. Let us consider the soft
operator,

$\displaystyle S^{abea^{\prime}b^{\prime}e^{\prime}}=\langle
0|(Y^{a^{\prime}b^{\prime}3d}\mathbf{v}\cdot\mathbfcal{B}_{s}^{e^{\prime}})^{\dagger}\mathcal{M}|X_{s}\rangle\langle
X_{s}|Y^{ab3d}\mathbf{v}\cdot\mathbfcal{B}_{s}^{e}|0\rangle,$ (115)

where $\mathcal{M}$ is the measurement performed on the soft operator and
$X_{s}$ are the soft modes. This soft operator now has 6 free indices which
must be contracted into the DM wavefunction factor. Since the soft final state
is completely inclusive, we can simplify the operator as

$\displaystyle S^{abea^{\prime}b^{\prime}e^{\prime}}=|\mathbf{v}|^{2}\langle
0|(Y^{a^{\prime}b^{\prime}3d}\mathbfcal{B}_{s}^{e^{\prime}i})^{\dagger}\mathcal{M}|X_{s}\rangle\langle
X_{s}|Y^{ab3d}\mathbfcal{B}_{s}^{ei}|0\rangle.$ (116)

We can now move the $|\mathbf{v}|^{2}$ factor to the DM wavefunction and
instead redefine a soft operator

$\displaystyle S^{abea^{\prime}b^{\prime}e^{\prime}}=\langle
0|(Y^{a^{\prime}b^{\prime}3d}\mathbfcal{B}_{s}^{e^{\prime}i})^{\dagger}\mathcal{M}|X_{s}\rangle\langle
X_{s}|Y^{ab3d}\mathbfcal{B}_{s}^{ei}|0\rangle.$ (117)

The only term that contributes at one loop is the $\mathbfcal{B}$ since it is
0 at ${\cal O}(\alpha_{\scriptscriptstyle W}^{0})$. So we can set all other
Wilson lines to their tree level values. Explicitly,

$\displaystyle
Y^{ab3d}\big{|}_{\textrm{tree}}=(Y_{v}^{fa}Y_{n}^{f3})(Y_{v}^{gb}Y_{\bar{n}}^{gd})\big{|}_{\textrm{tree}}=\delta^{a3}\delta^{bd}.$
(118)

Thus, our operator at one loop becomes

$\displaystyle S^{abea^{\prime}b^{\prime}e^{\prime}}_{1\mathchar
45\relax\mathrm{loop}}=\delta^{a3}\delta^{a^{\prime}3}\delta^{bb^{\prime}}\langle
0|(\mathbfcal{B}_{s}^{e^{\prime},i})^{\dagger}\mathcal{M}|X_{s}\rangle\langle
X_{s}|\mathbfcal{B}_{s}^{e,i}|0\rangle.$ (119)

We will only calculate this operator to one loop to elucidate its properties.
The one-loop integrand up to an overall factor take the form,

$\displaystyle I=$ $\displaystyle 2\delta^{ee^{\prime}}g_{\scriptscriptstyle
W}^{2}\int\frac{d^{d}k}{(2\pi)^{(d-1)}}\frac{\delta^{+}(k^{2}-M_{\chi}^{2})\delta(q^{+}-k^{+})}{k_{0}^{2}}$
(120) $\displaystyle-$ $\displaystyle\delta^{ee^{\prime}}g_{\scriptscriptstyle
W}^{2}\int\frac{d^{d}k}{(2\pi)^{(d-1)}}\frac{\delta^{+}(k^{2}-M_{\chi}^{2})\delta(q^{+}-k^{+})k^{2}}{k_{0}^{4}}.$

where $q^{+}$ is the contribution to the final photon momentum from the soft
function. The second term is proportional to $M_{\chi}^{2}$ and gives a power
correction in $M_{\chi}^{2}/(q^{+})^{2}$ and hence can be ignored. The first
term does not give a UV divergence, but will contribute a log which will be
relevant for stage 2 of the EFT. In detail,

$2\delta^{ee^{\prime}}g_{\scriptscriptstyle
W}^{2}\int\frac{d^{d}k}{(2\pi)^{(d-1)}}\frac{\delta^{+}(k^{2}-M_{\chi}^{2})\delta(q^{+}-k^{+})}{k_{0}^{2}}=2\delta^{ee^{\prime}}\frac{\alpha_{\scriptscriptstyle
W}}{\pi}\frac{q^{+}}{(q^{+})^{2}+M_{\chi}^{2}}.$ (121)

Going to Laplace space and expanding out in the limit $M_{\chi}\rightarrow 0$,
we have

$\displaystyle I=-2\delta^{ee^{\prime}}\frac{\alpha_{\scriptscriptstyle
W}}{\pi}\ln(M_{\chi}se^{\gamma_{E}}).$ (122)

This result has an IR divergence. At first glance, this may not be that
surprising since all our soft operators in the direct-channel annihilation
also were IR divergent. In those cases, we could trace back the IR divergence
to the violation of the KLN theorem due to the semi-inclusive nature of the
final state. In this case however, the interesting point is that, even when
the final state is completely inclusive, i.e., we do not constrain the final
state to be just a photon, our soft function is still IR divergent. Here we
can trace this to the exclusive nature of the initial state where we demand
that our operator support an $L=1$, $S=0$ state. This forces the emission of
the 3rd gauge boson which has no virtual counterpart, leading to an IR
divergence. This is very similar to the IR divergence that appears in the
computation of PDFs in QCD.

##  Appendix C Unstable Particle Effective Theory

In this section, we justify our use of Eq. (61) for computing the decay rate
of bound states. Let us now look at the effective theory for resonances
systematically. In the literature, this is referred to as unstable particle
effective theory (for a review see Ref. Beneke:2015vfa ).

To begin with, if we have an intermediate resonance state we expect a
propagator in our amplitude of the form

$\displaystyle D=\frac{i}{p^{2}-M_{*}^{2}},$ (123)

where $M_{*}$ is a complex pole. If we write $M_{*}=M+i\Gamma/2$, then we have
the result

$\displaystyle D=\frac{i}{p^{2}-M^{2}-i\Gamma M+\Gamma^{2}/4}.$ (124)

We wish to work in a regime of narrow width, i.e. $p^{2}-M^{2}\sim\Gamma M\ll
M^{2}$, so that the propagator is well approximated by

$\displaystyle D\simeq\frac{i}{p^{2}-M^{2}-i\Gamma M},$ (125)

and we want to develop an effective theory with an expansion in the small
parameter $\lambda=\Gamma/M$. For the case of inclusive decays of our DM bound
state, $\Gamma M_{\chi}\sim\alpha_{\scriptscriptstyle W}^{5}M_{\chi}^{2}$; we
are interested in the inclusive decay rate, so there are no large logarithms
and the perturbative cross section begins at $\alpha_{\scriptscriptstyle
W}^{2}$, and due to the nontrivial wavefunction, we have an additional factor
of at least $\alpha_{\scriptscriptstyle W}^{3}$. The effective coupling thus
scales as $\lambda\sim\alpha_{\scriptscriptstyle W}^{5}\ll 1$.

The hard scale in this process is just the resonance mass, which in our case
is simply $\sim$$M_{\chi}$. We can now treat this as an HQET theory writing
$p=M_{\chi}v+k$, where $v$ is the four-velocity of the resonance with
$v^{2}=1$ and $k$ is the residual momentum. Given the scaling above of
$p^{2}-M^{2}\sim\Gamma M$, we can immediately see that $k\sim\Gamma$, so that
we have a soft mode $k^{\mu}\sim(\Gamma,\Gamma,\Gamma)\equiv
M_{\chi}(\alpha_{\scriptscriptstyle W}^{5},\alpha_{\scriptscriptstyle
W}^{5},\alpha_{\scriptscriptstyle W}^{5})$. Obviously, the question is how
does this scale relate to the mass scale $m_{\scriptscriptstyle W}\sim
vM_{\chi}$ that we already have. Now the $\alpha_{\scriptscriptstyle W}$ in
the decay rate is evaluated at the scale $M_{\chi}$ while the
$\alpha_{\scriptscriptstyle W}$ in the HQET scaling is at the scale
$m_{\scriptscriptstyle W}$; however, we can see by the equations relating
$\alpha_{\scriptscriptstyle W}$ at the two scales that they are parametrically
of the same order. If that is the case, then our mode has $k^{\mu}\sim
M_{\chi}(\alpha_{\scriptscriptstyle W}^{5},\alpha_{\scriptscriptstyle
W}^{5},\alpha_{\scriptscriptstyle W}^{5})$ with $k^{2}\ll
m_{\scriptscriptstyle W}^{2}$ and hence can only be populated by a massless
mode such as the photon.

Given the fact that in our case $\lambda\sim\alpha_{\scriptscriptstyle
W}^{5}$, any corrections of the order $\lambda^{1}$ are minute and will be
sub-dominant in any error band at the accuracy we are aiming for. So we can
safely work at leading order in $\lambda$. Following Beneke:2015vfa , it is
clear that at this order the only term that exists is the propagator with the
1PI self-energy corrections. Any communication between the production and
decay states via radiative emissions of our mode $k^{\mu}$ only occurs at
$\mathcal{O}(\lambda)$ and hence is severely suppressed. Additionally, since
$\Gamma\sim\alpha_{W}^{5}M_{\chi}$, but the splitting between bound states,
$\Delta E_{n}\sim\alpha_{W}^{2}M_{\chi}$, any interference between bound
states is also subleading. Therefore, we can safely ignore any radiative
corrections by this mode. This suppression is a manifestation of the length
separation in space-time between the process of production and decay.

Now let us look at a cross-section for production of a resonance and its
subsequent decay. We assume that we are in a regime $p^{2}-M^{2}\sim\Gamma M$,
where $p$ is the momentum of the intermediate resonance state and $M$ is the
real part of the pole. Let $N$ be our initial scattering state that will
create the resonance, and we focus on the cross section to produce an observed
final state $f$ and an ultrasoft photon $\gamma_{\textrm{us}}$ indicating that
a bound state was formed. In detail, the differential cross section is

$\displaystyle\frac{d\sigma}{dz}=\frac{1}{\mathcal{N}}\int
d\Pi_{\gamma}d\Pi_{f}|{\cal M}(N\rightarrow
f+\gamma_{\textrm{us}})|^{2}\delta^{(4)}(p_{\gamma}+p_{f}-p_{N})\mathcal{M}_{z}(f),$
(126)

where ${\cal M}_{z}$ is the measurement (in this case the photon energy)
function on the final state particles, and ${\cal N}$ is a normalizing
kinematic factor. Since the photon emitted during the bound state formation is
an ultrasoft photon, there is no measurement on it (via a multipole expansion
of the measurement function) and its phase space is integrated over fully. The
amplitude squared will contain the squared propagator for the resonance,

$\displaystyle J=\frac{1}{(p^{2}-M^{2})^{2}+\Gamma^{2}M^{2}}.$ (127)

Now, the key point is that if we are not interested in the details of the
variation of the cross section near resonance, and we are sufficiently
inclusive over $p^{2}$ around the resonance (by a value $\gg\Gamma M$), then
we can make the following substitution

$\displaystyle\lim_{\Gamma/M\rightarrow 0}J\rightarrow\frac{\pi}{\Gamma
M}\delta(p^{2}-M^{2}).$ (128)

This substitution is true only in the distribution sense, i.e. under the
integral which at least encompasses the region of the size of the width about
the resonance. This is the narrow width approximation. This substitution then
puts the intermediate resonance on-shell. The cross section can then be
written as

$\displaystyle\frac{d\sigma}{dz}$
$\displaystyle=\frac{1}{\mathcal{N}}\frac{\pi}{\Gamma M}\int
d\Pi_{\gamma}d\Pi_{f}|{\cal M}(N\rightarrow
B(p)+\gamma_{\textrm{us}})|^{2}|{\cal M}(B(p)\rightarrow f)|^{2}$ (129)
$\displaystyle\times\delta^{(4)}(p_{\gamma}+p_{f}-p_{N})\delta(p^{2}-M^{2})\mathcal{M}_{z}(f).$

Here $B(p)$ represents a bound state with momentum $p$, and again this result
is true at leading order in $\Gamma/M$, which forbids any communication
between the production and decay states. If we then insert a factor of unity,
$1=\int d^{4}p\,\delta^{(4)}(p-p_{f})$, the result can be rearranged to yield,

$\displaystyle\frac{d\sigma}{dz}=$
$\displaystyle\frac{1}{\mathcal{N}}\frac{\pi}{M}\Big{[}\int
d\Pi_{\gamma}d\Pi_{R}|{\cal M}(N\rightarrow
B(p)+\gamma_{\textrm{us}})|^{2}\delta^{(4)}(p_{\gamma}+p-p_{N})\Big{]}$ (130)
$\displaystyle\times$ $\displaystyle\frac{1}{\Gamma}\Big{[}\int d\Pi_{f}|{\cal
M}(B(p)\rightarrow f)|^{2}\delta^{(4)}(p-p_{f})\mathcal{M}_{z}(f)\Big{]}$
$\displaystyle=$ $\displaystyle\sigma(N\rightarrow
B+\gamma_{\textrm{us}})\frac{1}{\Gamma}\frac{d\Gamma_{B\rightarrow f}}{dz},$

which is simply the product of the production cross section and the
differential branching ratio. The above separation holds where the process
proceeds solely through the long lived bound state, but in practice the result
should be summed over all bound states in the spectrum, as well as the direct
annihilation case where no bound states are formed. Applied to our specific
case we arrive at Eq. (61).

##  Appendix D Proof of the Wilson Line Identity

In this section, we prove the identity involving soft Wilson lines used in Eq.
(10) that eventually leads to a universal factorization of the IR physics in
terms of soft and jet functions independent of the representation. The
property that we wish to show is

$\displaystyle S_{v}T^{a}S_{v}^{\dagger}=T^{a^{\prime}}S_{v}^{a^{\prime}a},$
(131)

where $T$ is a generator in an arbitrary representation (we used $\mathbf{5}$
for the quintuplet), and on the left hand side we have two Wilson lines in the
same representation, whereas on the right it is in the adjoint. (In the main
text we used $Y_{v}$ for the latter, we keep all as $S$ here for notational
convenience.) In the main text we actually used
$S_{v}^{\dagger}T^{a}S_{v}=S_{v}^{aa^{\prime}}T^{a^{\prime}}$, although this
follows from the above by applying various inverses. In position space, the
Wilson lines are defined as

$\displaystyle S_{v}(x)$ $\displaystyle=Pe^{ig\int_{-\infty}^{x}ds\,v\cdot
A_{s}(vs)}=Pe^{ig\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot y\,v\cdot
A_{s}(\bar{v}\cdot y)},$ (132) $\displaystyle S_{v}^{\dagger}(x)$
$\displaystyle=\bar{P}e^{-ig\int_{-\infty}^{x}ds\,v\cdot A_{s}(vs)}.$

where $v\cdot A_{s}=v\cdot A_{s}^{a}T^{a}$, with $T^{a}$ in the appropriate
representation for the Wilson line, whilst $P$ is path ordering and $\bar{P}$
indicates anti-path ordering. The variable $s$ parametrizes the path along the
light cone direction $n$ from $x$ to $-\infty$. The statement in Eq. (131) is
a generalization of the identity applied for QCD Becher:2014oda  for other
group representations. The soft Wilson line obeys the equation

$\displaystyle\frac{d}{d\bar{v}\cdot x}S_{v}(x)$ $\displaystyle=ig\,v\cdot
A_{s}(\bar{v}\cdot x)S_{v}(x),$ (133) $\displaystyle\frac{d}{d\bar{v}\cdot
x}S^{\dagger}_{v}(x)$ $\displaystyle=-S_{v}^{\dagger}(x)ig\,v\cdot
A_{s}(\bar{v}\cdot x).$

Coming back to our question, let us define
$U^{a}(x)=S_{v}(x)T^{a}S_{v}^{\dagger}(x)$. We can then immediately see

$\displaystyle\frac{d}{d\bar{v}\cdot x}U^{a}(x)=\left[igv\cdot
A_{s}(\bar{v}\cdot x),U^{a}(x)\right]\\!.$ (134)

We will solve this equation by recursion, order by order in the coupling $g$
to build up the full solution. The tree level result is simply
$U^{a(0)}(x)=T^{a}=T^{a^{\prime}}\delta^{a^{\prime}a}$. At the next order,

$\displaystyle U^{a(1)}(x)$ $\displaystyle=\int_{-\infty}^{\bar{v}\cdot
x}d\bar{v}\cdot y\left[ign\cdot A_{s}(\bar{v}\cdot y),T^{a}\right]$ (135)
$\displaystyle=ig\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot y\,v\cdot
A_{s}^{b}(\bar{v}\cdot y)if^{baa^{\prime}}T^{a^{\prime}}$
$\displaystyle=T^{a^{\prime}}ig\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot
y\,v\cdot(A_{s}(\bar{v}\cdot y))^{a^{\prime}a}.$

At ${\cal O}(g^{2})$,

$\displaystyle U^{a(2)}(x)$ $\displaystyle=\int_{-\infty}^{\bar{v}\cdot
x}d\bar{v}\cdot y_{1}\int_{-\infty}^{\bar{v}\cdot y_{1}}d\bar{v}\cdot
y_{2}\left[igv\cdot A_{s}(\bar{v}\cdot y_{1}),\left[igv\cdot
A_{s}(\bar{v}\cdot y_{2}),T^{a}\right]\right]$ (136)
$\displaystyle=\frac{1}{2}\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot
y_{1}\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot
y_{2}P\Big{\\{}\left[igv\cdot A_{s}(\bar{v}\cdot y_{1}),\left[igv\cdot
A_{s}(\bar{v}\cdot y_{2}),T^{a}\right]\right]\Big{\\}}$
$\displaystyle=T^{a^{\prime}}\frac{(ig)^{2}}{2}\int_{-\infty}^{\bar{v}\cdot
x}d\bar{v}\cdot y_{1}\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot
y_{2}P\Big{\\{}\left(v\cdot A_{s}(\bar{v}\cdot y_{1})v\cdot A_{s}(\bar{v}\cdot
y_{2})\right)^{a^{\prime}a}\Big{\\}}.$

From here, we can see that the $n^{th}$ term will be

$\displaystyle U^{a(n)}(x)$
$\displaystyle=\frac{1}{n!}\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot
y_{1}\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot
y_{2}\ldots\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot y_{n}$ (137)
$\displaystyle\times P\Big{\\{}\left[igv\cdot A_{s}(\bar{v}\cdot
y_{1}),\left[igv\cdot A_{s}(\bar{v}\cdot y_{2}),\left[\ldots\left[igv\cdot
A_{s}(\bar{v}\cdot y_{n}),T^{a}\right]\right]\right]\right]\Big{\\}}$
$\displaystyle=T^{a^{\prime}}\frac{(ig)^{n}}{n!}\int_{-\infty}^{\bar{v}\cdot
x}d\bar{v}\cdot y_{1}\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot
y_{2}\ldots\int_{-\infty}^{\bar{v}\cdot x}d\bar{v}\cdot y_{n}$
$\displaystyle\times P\Big{\\{}\left(v\cdot A_{s}(\bar{v}\cdot y_{1})v\cdot
A_{s}(\bar{v}\cdot y_{2})\ldots v\cdot A_{s}(\bar{v}\cdot
y_{n})\right)^{a^{\prime}a}\Big{\\}}.$

Summing to all orders then proves our result.

##  Appendix E Subtle Signs in the Bound-state Formation and Decay
Calculations

In the bulk of the paper, we have freely used results from Ref. Harz:2018csl
that are written in terms of two-body states of the form $|ij\rangle$.
However, in general our 2-body states for non-identical particles will in fact
be combinations of the form
$\frac{1}{\sqrt{2}}(|ij\rangle+(-1)^{L+S}|ji\rangle)$. (This convention choice
is also discussed in the context of Sommerfeld enhancement in Ref.
Beneke:2014gja , where these two approaches are labeled “method-1” and
“method-2”; we largely adopt “method-2”, where we treat $|ij\rangle$ and
$|ji\rangle$ as components of a single state, rather than tracking them
separately.) The factor of $(-1)^{L+S}$ arises from a factor of $(-1)^{S+1}$
from the behavior of the spin configuration under particle exchange, a factor
of $(-1)^{L}$ from the parity of the spatial wavefunction, and a factor of
$(-1)$ from the exchange of two fermions.

This means that when considering a transition of the form
$|ij\rangle\rightarrow|i^{\prime}j^{\prime}\rangle$, we also need to include
transitions between the $|ji\rangle$ and $|j^{\prime}i^{\prime}\rangle$
states. In many cases this does not make a difference and it is adequate to
represent states purely by one component $|ij\rangle$. For example, if the
$|ij\rangle\rightarrow|i^{\prime}j^{\prime}\rangle$ and
$|ji\rangle\rightarrow|j^{\prime}i^{\prime}\rangle$ processes have equal
rates, but $|ij\rangle\rightarrow|j^{\prime}i^{\prime}\rangle$ and
$|ji\rangle\rightarrow|i^{\prime}j^{\prime}\rangle$ are forbidden (for
example, this occurs if $|ij\rangle=|0,++\rangle$), then the combined rate is
the same as what one would obtain from purely considering the
$|ij\rangle\rightarrow|i^{\prime}j^{\prime}\rangle$ process. However, this
behavior is not universal.

As an example of a case where this matters, consider the transition between
$Q=1$ bound state components $|\\!+0\rangle\rightarrow|\\!+0\rangle$. Writing
out the individual components of these states (labeled as $23$ and $32$,
following the notation of Sec. 3),262626One might be tempted to write the 23
state as $|\\!+0\rangle$ and 32 as $|0+\rangle$. However, we reserve
$|\\!+0\rangle$ for the full quantum state such that
$|\\!+0\rangle=(|2\,3\rangle+(-1)^{L+S}|3\,2\rangle)/\sqrt{2}$. the full
matrix element should be:

$\displaystyle\mathcal{M}=\frac{1}{2}$
$\displaystyle\left(\mathcal{M}_{22,33}+(-1)^{(L+S)_{i}}\mathcal{M}_{32,23}+(-1)^{(L+S)_{f}}\mathcal{M}_{23,32}\right.$
(138)
$\displaystyle\left.+(-1)^{(L+S)_{i}+(L+S)_{f}}\mathcal{M}_{33,22}\right)\\!.$

Now we can write $(L+S)_{f}=(L+S)_{i}+1\,\textrm{(mod 2)}$ for dipole
transitions, and consequently:

$\mathcal{M}=\frac{1}{2}\left(\mathcal{M}_{22,33}-\mathcal{M}_{33,22}+(-1)^{(L+S)_{i}}\left(\mathcal{M}_{32,23}-\mathcal{M}_{23,32}\right)\right)\\!.$
(139)

Now as calculated in Sec. 3, if $\psi_{i}$ and $\psi_{f}$ denote the initial-
and final-state wavefunctions, we have:

$\displaystyle\mathcal{M}^{3}_{22,33}$
$\displaystyle=i\sqrt{2^{6}\pi\alpha_{\text{rad}}M_{\chi}}\left[(T^{3})_{22}-(T^{3})_{33}\right]\int
d^{3}\mathbf{r}\,\psi_{f}^{*}\nabla\psi_{i}$ (140)
$\displaystyle=i\sqrt{2^{6}\pi\alpha_{\text{rad}}M_{\chi}}\int
d^{3}\mathbf{r}\,\psi_{f}^{*}\nabla\psi_{i},{}$
$\displaystyle\mathcal{M}^{3}_{33,22}$
$\displaystyle=i\sqrt{2^{6}\pi\alpha_{\text{rad}}M_{\chi}}\left[(T^{3})_{33}-(T^{3})_{22}\right]\int
d^{3}\mathbf{r}\,\psi_{f}^{*}\nabla\psi_{i}{}$
$\displaystyle=i\sqrt{2^{6}\pi\alpha_{\text{rad}}M_{\chi}}\int
d^{3}\mathbf{r}\,(-\psi_{f}^{*}\nabla\psi_{i}),{}$
$\displaystyle\mathcal{M}^{3}_{23,32}$
$\displaystyle=i\sqrt{2^{6}\pi\alpha_{\text{rad}}M_{\chi}}\left\\{\left[-i((T^{1})_{32}(T^{2})_{23}-(T^{2})_{32}(T^{1})_{23}\right]M_{\chi}\alpha_{\text{NA}}\int
d^{3}\mathbf{r}\,\hat{\mathbf{r}}\psi_{f}^{*}\psi_{i}\right\\}{}$
$\displaystyle=i\sqrt{2^{6}\pi\alpha_{\text{rad}}M_{\chi}}\left\\{-3M_{\chi}\alpha_{\text{NA}}\int
d^{3}\mathbf{r}\,\hat{\mathbf{r}}\psi_{f}^{*}\psi_{i}\right\\}\\!,{}$
$\displaystyle\mathcal{M}^{3}_{32,23}$
$\displaystyle=i\sqrt{2^{6}\pi\alpha_{\text{rad}}M_{\chi}}\left\\{\left[-i((T^{1})_{23}(T^{2})_{32}-(T^{2})_{23}(T^{1})_{32}\right]M_{\chi}\alpha_{\text{NA}}\int
d^{3}\mathbf{r}\,\hat{\mathbf{r}}\psi_{f}^{*}\psi_{i}\right\\}{}$
$\displaystyle=i\sqrt{2^{6}\pi\alpha_{\text{rad}}M_{\chi}}\left\\{3M_{\chi}\alpha_{\text{NA}}\int
d^{3}\mathbf{r}\,\hat{\mathbf{r}}\psi_{f}^{*}\psi_{i}\right\\}\\!,{}$

recalling that the “3” superscript is for $\gamma,\,Z$ emission, depending on
the value of $\alpha_{\text{rad}}$, the coupling of the emitted boson to the
charged particle that radiated it. We also recall that the terms with
$\alpha_{\text{NA}}$ correspond to emission off the virtual particles sourcing
the potential. The $\alpha_{\text{NA}}$ factor is thus the coupling of the
virtual line to the WIMPs. For the case of capture by a radiated $\gamma$,
$\alpha_{\text{rad}}=\alpha_{\text{em}}$ and $\alpha_{\text{NA}}=\alpha_{W}$.

Thus, overall we have $\mathcal{M}_{32,23}=-\mathcal{M}_{23,32}$ and
$\mathcal{M}_{33,22}=-\mathcal{M}_{22,33}$, and consequently:

$\displaystyle\mathcal{M}$
$\displaystyle=\mathcal{M}_{22,33}+(-1)^{(L+S)_{i}}\mathcal{M}_{32,23}{}$
$\displaystyle=i\sqrt{2^{6}\pi\alpha_{\text{rad}}M_{\chi}}\left[\int
d^{3}\mathbf{r}\,\psi_{f}^{*}\nabla\psi_{i}+(-1)^{(L+S)_{i}}3M_{\chi}\alpha_{\text{NA}}\int
d^{3}\mathbf{r}\,\hat{\mathbf{r}}\psi_{f}^{*}\psi_{i}\right]\\!.$ (141)

The $(-1)^{(L+S)_{i}}$ factor is obtained by treating the different components
correctly, and is required to ensure the correct behavior of the matrix
element under time reversal.

##  Appendix F Analytic Approximate Results for Annihilation and Bound-state
Formation

In this final appendix we provide analytic estimates for the annihilation and
bound state formation rate of DM. We consider the quintuplet case of interest
first, followed by providing equivalent results for a general representation.

###  F.1 Results for the quintuplet

In the limit of unbroken SU(2), the wavefunctions and their integrals, and
hence the bound-state capture rate, can be computed analytically in the low-
velocity limit. In this regime the Sommerfeld enhancement can also be computed
analytically. These results can also be applied (approximately, and with
caveats we will discuss below) to the case where SU(2) is broken but the DM
mass is very heavy relative to the symmetry breaking scale. These calculations
can be useful both as a cross-check on our numerical results, and to develop
intuition for which channels are likely to dominate the overall annihilation
signal. As such, we present the details of these analytic calculations below,
beginning with the DM in the quintuplet representation.

####  F.1.1 Capture and annihilation rates

As an opening example, let us estimate the spin-averaged capture rate into the
spin-triplet ground state via photon emission. The total cross section for
this process is given by (see App. C of Ref. Asadi:2016ybp ),272727This
expression as written includes capture only from the components of the
incoming state that experience an attractive potential; we expect the
contribution from repulsed incoming states to be suppressed, due to their
small overlap with the bound states.

$\displaystyle\sigma v$
$\displaystyle=\frac{3}{2}\times\frac{2^{8}\pi\alpha\,k}{3}\left|\sum_{i}({\bf
I}\cdot\eta_{i})\alpha_{\scriptscriptstyle W}(\alpha_{\scriptscriptstyle
W}\lambda_{f}M_{\chi}/2)^{-3/2}e^{-2\lambda_{i}/\lambda_{f}}e^{\pi\alpha_{\scaleto{W}{3.5pt}}\lambda_{i}/(2v)}\Gamma(1-i\alpha_{\scriptscriptstyle
W}\lambda_{i}/v)\right.{}$
$\displaystyle\left.\times\eta_{f}^{\dagger}\left[\lambda_{i}\hat{C}_{1}+\hat{C}_{2}\lambda_{i}/\lambda_{f}\right]\eta_{i}\right|^{2}\\!,$
(142)

The notation we employ follows Ref. Asadi:2016ybp ; $\eta_{i}$ and $\eta_{f}$
are potential eigenvectors which for the quintuplet are given by (in our
basis) $\eta_{f}=\\{-2,1,0\\}/\sqrt{5}$, $i=1,2$ with
$\eta_{1}=\\{\sqrt{2},-\sqrt{2},1\\}/\sqrt{5}$,
$\eta_{2}=\\{-2,-1,\sqrt{2}\\}/\sqrt{7}$, with corresponding attractive
eigenvalues $\lambda_{f}=5$, $\lambda_{1}=6$, $\lambda_{2}=3$. The $I$ vector
describes the fraction of the incoming plane wave in each state; we will
choose $I=\\{0,0,1\\}$, as the state asymptotes to two noninteracting, neutral
DM particles. The energy of the outgoing photon is $k$, which at low
velocities can be approximated as the binding energy of the ground state,
$\lambda_{f}^{2}\alpha_{\scriptscriptstyle W}^{2}M_{\chi}/4$. Lastly, the
$\hat{C}_{1}$ and $\hat{C}_{2}$ matrices describe the couplings between the
different components of the initial and final states; for capture via photon
(or $Z$) emission, they take the form:

$\hat{C}_{1}=\begin{pmatrix}2&0&0\\\ 0&1&0\\\
0&0&0\end{pmatrix}\\!,\quad\hat{C}_{2}=\begin{pmatrix}0&2&0\\\
-2&0&3\sqrt{2}\\\ 0&-3\sqrt{2}&0\end{pmatrix}\\!.$ (143)

In the wino and positronium cases studied in Ref. Asadi:2016ybp , there was
only one eigenstate that experienced an attractive initial-state potential,
and so the sum in Eq. (142) was trivial. There is a simple expression for
$|e^{\pi\alpha_{\scaleto{W}{3.5pt}}\lambda_{i}/(2v)}\Gamma(1-i\alpha_{\scriptscriptstyle
W}\lambda_{i}/v)|^{2}\simeq 2\pi\alpha_{\scriptscriptstyle W}\lambda_{i}/v$ in
the limit of small $v$ (for positive $\lambda_{i}$), which scales purely as
$1/v$, and consequently in those cases (wino and positronium) $\sigma v$ had a
simple $1/v$ scaling at low relative velocities. However, in the quintuplet
case, we see there are multiple terms in the sum that can interfere with each
other, and so even in the limit of unbroken SU(2), we expect there to be a
non-trivial velocity dependence in the capture cross section.

Similarly, the Sommerfeld factors can be read off from the components of the
scattering-state wavefunction at the origin, and in the unbroken limit this
wavefunction has the form
$\sum_{i}(I\cdot\eta_{i})\eta_{i}\,\phi(\lambda_{i}\alpha_{\scriptscriptstyle
W},r)$, where $\phi(\alpha,r)$ is the solution to the scalar Schrodinger
equation with an attractive Coulomb potential with coupling $\alpha$. In
principle this sum runs over both positive and negative eigenvalues of the
potential (corresponding to both attracted and repulsed eigenstates), but for
low velocities we expect the contribution of the eigenstates experiencing a
repulsive interaction to be very small. Nonetheless, where (as in the
quintuplet case) there are two eigenstates experiencing an attractive
potential, the value of the wavefunction at the origin (and hence the
Sommerfeld factors) will experience a non-trivial interference between the two
contributions. This can give rise to a velocity dependence differing from the
simpler case where there is only one attractive eigenstate. Similar
interference effects can be seen in the form of rapid changes in spectrum with
respect to $M_{\chi}$ in the case of broken SU(2) symmetry, where the
interference occurs between the various Sommerfeld factors with resonances at
different positions (as discussed in the context of Fig. 9). In contrast, the
manifestation of the eigenstate interference identified here persists in the
SU(2)-symmetric limit and does not require any resonance structure, only
differing (velocity-dependent) phases between the interfering contributions.

However, as noted in Ref. Schutz:2014nka ; Asadi:2016ybp , at low velocities
the system is often in an “adiabatic” regime where the incoming particle
wavefunction evolves such that at short distances it has complete overlap with
the eigenvector with the largest-magnitude attractive eigenvalue. The
criterion for this behavior is roughly $v\lesssim\delta/m_{\scriptscriptstyle
W}$, where $\delta$ is the mass splitting between the states; for $\delta=164$
MeV, we expect this behavior to hold roughly for $v\lesssim 2\times 10^{-3}$,
i.e. for Milky-Way-scale velocities and lower. Note that this criterion is
independent of the DM mass, so even when the DM is very heavy and the ratio
$m_{\scriptscriptstyle W}/M_{\chi}$ is small, the effect of SU(2) breaking can
still be seen in the presence of this adiabatic regime. In this case, the
interference will be suppressed for both bound-state formation and Sommerfeld
enhancement, with only the $i=1$, $\lambda_{i}=6$ term contributing
significantly, and with the coefficient $I\cdot\eta_{i}$ replaced with
$\delta_{i1}$.

This is an important simplifying approximation within its regime of validity.
Note that the presence of this regime relies on SU(2) being broken, and also
on low velocity ($v\lesssim\delta/m_{\scriptscriptstyle W}$); it will not
appear if an unbroken symmetry ensures the degeneracy of the mass eigenstates,
and it will also not generally be relevant in the early universe (e.g. for
relic density calculations) where velocities are much higher. However, it is
well-suited to the case of indirect detection in the Milky Way halo.

For example, within this approximation, we obtain the spin-averaged capture
rate to the ground state as:

$\sigma v\simeq\frac{2^{8}}{5}\frac{\pi\alpha\alpha_{\scriptscriptstyle
W}}{M_{\chi}^{2}}\times\frac{3^{3}\cdot
2^{9}}{5^{2}}e^{-24/5}\frac{\pi\alpha_{\scriptscriptstyle
W}}{v}=\frac{2^{17}\cdot
3^{3}}{5^{3}}e^{-24/5}\frac{\pi^{2}\alpha\alpha_{\scriptscriptstyle
W}^{2}}{M_{\chi}^{2}v}\simeq\frac{233\pi^{2}}{v}\frac{\alpha\alpha_{\scriptscriptstyle
W}^{2}}{M_{\chi}^{2}}.$ (144)

Here we have employed the low-velocity approximation
$|e^{\pi\alpha_{\scaleto{W}{3.5pt}}\lambda_{i}/(2v)}\Gamma(1-i\alpha_{\scriptscriptstyle
W}\lambda_{i}/v)|^{2}\simeq 2\pi\alpha_{\scriptscriptstyle W}\lambda_{i}/v$.

In the same regime, where the initial state rotates into the most-attracted
eigenstate, the $s$-wave direct annihilation cross section to gauge bosons can
be computed as,

$\sigma v\simeq\frac{720\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}.$ (145)

In this unbroken limit, the effective branching ratio to the line (i.e. to
$\gamma\gamma$ \+ half the branching ratio to $\gamma Z$) should be given by
$(s_{\scriptscriptstyle W}^{4}+s_{\scriptscriptstyle
W}^{2}c_{\scriptscriptstyle W}^{2})/3=s_{\scriptscriptstyle W}^{2}/3$, so the
line cross section should be:

$(\sigma v)_{\text{line}}\simeq\frac{240\pi^{2}\alpha_{\scriptscriptstyle
W}^{2}\alpha}{M_{\chi}^{2}v}.$ (146)

As for the bound-state formation, at higher velocities ($v\gtrsim 2\times
10^{-3}$), we expect to see the onset of interference between the
contributions from the two attracted eigenstates, resulting in a non-monotonic
dependence of the cross section on velocity even in the $s$-wave case. This
behavior, and its onset at roughly Milky Way-scale velocities, can be observed
in Fig. 10. (Note that the velocity dependence can also be suppressed if the
DM mass is small enough that the Sommerfeld enhancement is fully saturated,
such that the velocity dependence of the individual Sommerfeld factors is very
different from the case of unbroken SU(2) symmetry.)

This would suggest the cross-sections for capture (to the spin-triplet ground
state) and for annihilation producing a line should be very similar, at least
when this adiabatic approximation holds (numerical calculations indicate the
non-adiabatic cross section for bound-state capture, as estimated in Eq. 142,
can range between larger than the adiabatic result by a factor of $\sim$$2$
and smaller by a factor of $\sim$$5$, as $v$ is varied). However, for $v\ll
m_{\scriptscriptstyle W}/M_{\chi}$, we have the SU(2)-breaking effect that
$p$-wave processes should be parametrically suppressed by a factor of order
$(vM_{\chi}/m_{\scriptscriptstyle W})^{2}$, which at a 13.6 TeV mass can
suppress the $p$-wave capture cross section by $\sim$$2$ orders of magnitude.

We can also study the cross-section for spin-singlet $s\rightarrow p$ capture
to an $n=2$, $l=1$ state. In this case we have
$k=(25/16)\alpha_{\scriptscriptstyle W}^{2}M_{\chi}$, and

$\displaystyle\sigma v$ $\displaystyle=\frac{2^{13}\pi\alpha
k}{3^{3}}\frac{1}{\alpha_{\scriptscriptstyle
W}}M_{\chi}^{-3}\frac{1}{\lambda_{f}^{3}}\left|\sum_{i}({\bf
I}\cdot\eta_{i})e^{-4\lambda_{i}/\lambda_{f}}e^{\pi\alpha_{\scaleto{W}{3.5pt}}\lambda_{i}/(2v)}\Gamma(1-i\alpha_{\scriptscriptstyle
W}\lambda_{i}/v)\eta_{f}^{\dagger}\right.$ (147)
$\displaystyle\left.\times\left[\lambda_{i}\left(\frac{4\lambda_{i}}{\lambda_{f}}-3\right)\hat{C}_{1}+\hat{C}_{2}\left(3-12\frac{\lambda_{i}}{\lambda_{f}}+8\frac{\lambda_{i}^{2}}{\lambda_{f}^{2}}\right)\right]\eta_{i}\right|^{2}$
$\displaystyle=\frac{2^{9}\pi\alpha\alpha_{\scriptscriptstyle W}}{5\cdot
3^{3}M_{\chi}^{2}}\Bigg{|}\frac{3}{25}\sqrt{\frac{2}{5}}e^{-24/5}e^{3\pi\alpha_{\scaleto{W}{3.5pt}}/(2v)}\left(37e^{12/5}\Gamma(1-3i\alpha_{\scriptscriptstyle
W}/v)\right.$ $\displaystyle\hskip
184.9429pt\left.+89e^{3\pi\alpha_{\scaleto{W}{3.5pt}}/(2v)}\Gamma(1-6i\alpha_{\scriptscriptstyle
W}/v)\right)\Bigg{|}^{2},$

or in the adiabatic regime,

$\displaystyle\sigma v$
$\displaystyle\rightarrow\frac{2^{9}\pi\alpha\alpha_{\scriptscriptstyle
W}}{5\cdot
3^{3}M_{\chi}^{2}}\left|\frac{267}{25}\sqrt{2}e^{-24/5}e^{3\pi\alpha_{\scaleto{W}{3.5pt}}/(2v)}\Gamma(1-6i\alpha_{\scriptscriptstyle
W}/v)\right|^{2}$ (148) $\displaystyle=\frac{2^{12}\cdot
89^{2}\pi^{2}}{5^{5}\,v}e^{-48/5}\frac{\alpha\alpha_{\scriptscriptstyle
W}^{2}}{M_{\chi}^{2}}$
$\displaystyle\simeq\frac{0.70\pi^{2}}{v}\frac{\alpha\alpha_{\scriptscriptstyle
W}^{2}}{M_{\chi}^{2}}.$

We see that the scale for this cross-section is naturally two orders of
magnitude smaller than the previous ones, which arises from the various
numerical prefactors, primarily the factor of $e^{-48/5}$ compared to
$e^{-24/5}$ for the capture to the $n=1$ state, corresponding to a suppression
factor of $8\times 10^{-3}$. If these exponential terms were removed, the
other prefactors would differ by less than a factor of 3. (Numerical
calculations indicate the non-adiabatic cross section is larger than the
adiabatic one in this case, by factors between 1 and 3.6.) This is suggestive
that only the capture to the ground-state is likely to be comparable to direct
annihilation, and the capture from the $p$-wave initial-state component
suffers from a $v^{2}$ suppression once $v$ drops below $m_{\scriptscriptstyle
W}/M_{\chi}$, which renders it subdominant at the $\mathcal{O}(1\%)$ level for
our 13.6 TeV benchmark point.

This suppression for higher-$n$ capture also suggests that the contribution to
the endpoint hard photon spectrum from bound state formation and decay will be
suppressed, as these contributions are dominated by capture into states with
odd $L$ (and thus $n>1$) and $S=0$ that decay to $L=S=0$ states before
annihilating (see Sec. 4 for a more in-depth discussion). Capture to the
ground-state via emission of a dipole gauge boson changes $L$ by 1, thus
requiring an initial $L=1$ state (which must then have $S=1$ if it contains
identical DM particles), and $S=1$ states do not produce a leading-power
contribution to the endpoint spectrum when they decay.

For the $Q=1$ sector, let us again consider capture from the spin-triplet
$p$-wave incoming wave to the spin-triplet $s$-wave state. In the unbroken
limit the potential matrix for the final state takes the form:

$V=\begin{pmatrix}-2&\sqrt{6}\\\ \sqrt{6}&-3\end{pmatrix}\\!.$ (149)

Here the first row refers to the $++-$ state and the second to the $+\,0$
state. The transition matrices are now,

$\hat{C}_{1}=\begin{pmatrix}-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\\
0&-\frac{\sqrt{3}}{2}&\sqrt{\frac{3}{2}}\end{pmatrix}\\!,\quad\hat{C}_{2}=\begin{pmatrix}2\sqrt{2}&\sqrt{2}&0\\\
0&\sqrt{3}&0\end{pmatrix}\\!.$ (150)

In this case we must also replace the $\alpha$ prefactor in the cross section
with $\alpha_{\scriptscriptstyle W}$. The attractive eigenvalue for the final
state has $\lambda_{f}=5$, $\eta_{f}=\\{-\sqrt{2},\sqrt{3}\\}/\sqrt{5}$.
Again, the binding energy of the ground state is $k=\alpha_{\scriptscriptstyle
W}^{2}\lambda_{f}^{2}M_{\chi}/4=(25/4)\alpha_{\scriptscriptstyle
W}^{2}M_{\chi}$ (and since in the unbroken limit $m_{\scriptscriptstyle W}=0$,
we do not need to include a kinematic suppression for the $W$ mass). Then we
obtain for the unbroken limit:

$\displaystyle\sigma v=\frac{2^{7}3^{2}\pi\alpha_{\scriptscriptstyle
W}^{2}}{5^{4}M_{\chi}^{2}}$
$\displaystyle\left|16e^{-12/5}e^{6\pi\alpha_{\scaleto{W}{3.5pt}}/(2v)}\Gamma(1-6i\alpha_{\scriptscriptstyle
W}/v)\right.$ (151)
$\displaystyle\left.+7e^{-6/5}e^{3\pi\alpha_{\scaleto{W}{3.5pt}}/(2v)}\Gamma(1-3i\alpha_{\scriptscriptstyle
W}/v)\right|^{2}\\!,$

When we assume the adiabatic regime, the result reduces to,

$\sigma v=\frac{2^{17}\cdot
3^{3}\pi^{2}}{5^{3}}e^{-24/5}\frac{\alpha_{\scriptscriptstyle
W}^{3}}{vM_{\chi}^{2}}\simeq 233\frac{\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{vM_{\chi}^{2}}.$ (152)

This is for capture to the $Q=+1$ state; there is an equal rate for capture to
the $Q=-1$ state. Note the $\alpha_{\scriptscriptstyle W}$ prefactor (rather
than $\alpha$); including formation of the $Q=0$ state through $Z$ emission
(as well as photon emission) would similarly promote that capture rate to have
a prefactor of $\alpha_{\scriptscriptstyle W}$ rather than $\alpha$, for an
overall capture rate (summing across all three channels) of $\sim
700\pi^{2}\alpha_{\scriptscriptstyle W}^{3}/(M_{\chi}^{2}v)$, similar to the
full $s$-wave direct annihilation rate. The primary difference between this
$p\rightarrow s$ capture rate and the inclusive direct annihilation rate will
arise from velocity suppression of the $p\rightarrow s$ capture cross section
in the broken-SU(2) case (with this suppression being lifted in the truly
unbroken limit).

####  F.1.2 Presence of metastable bound states

The unbroken-SU(2) limit is also helpful for studying the question of whether
there could be $L>0$ states in the spectrum whose decays to more deeply bound
states are highly suppressed, leading them to decay through annihilation to SM
particles with a substantial branching ratio. States which are degenerate in
the unbroken limit are likely to remain close in energy as we reduce
$M_{\chi}$, and therefore decays between them will be suppressed (although if
this is decisive in determining whether a state is metastable, a more careful
analysis will be required). The $L+S$-even potential for the quintuplet has
two attractive eigenvalues, $Z=6$ and $Z=3$, whereas the $L+S$-odd potential
has a single attractive eigenvalue, $Z=5$. Thus for spin-singlet states
($S=0$) we expect $L$-even bound states with energies
$E_{n}/\alpha_{\scriptscriptstyle W}^{2}M_{\chi}=-9/n^{2},-2.25/n^{2}$, and
$L$-odd bound states with energies $E_{n}/\alpha_{\scriptscriptstyle
W}^{2}M_{\chi}=-6.25/n^{2}$. For spin-triplet states ($S=1$) we expect $L$-odd
states with energies $E_{n}/\alpha_{\scriptscriptstyle
W}^{2}M_{\chi}=-9/n^{2},-2.25/n^{2}$, and $L$-even states with energies
$E_{n}/\alpha_{\scriptscriptstyle W}^{2}M_{\chi}=-6.25/n^{2}$.

We first consider the case of $L$-odd states. Given a spin-singlet $L$-odd
state with $n>L>0$ (binding energy $6.25/n^{2}$), there should always be a
more deeply bound state with $n^{\prime}=n$, $L^{\prime}=L-1$ (binding energy
$9/n^{2}$), which is accessible through a dipole transition. So in the spin-
singlet case and Coulomb limit there should be no metastable states with
$L>0$. For the spin-triplet, the case is slightly more complicated, as given
an $L$-odd state with $n>L>0$ (binding energy $9/n^{2}$ or $2.25/n^{2}$), the
accompanying state with $n^{\prime}=n$, $L^{\prime}=L-1$ has binding energy
$6.25/n^{2}$. We see that the $L$-odd states with binding energies
$2.25/n^{2}$ will always have an accompanying more-deeply-bound $L$-even spin-
triplet state, to which they can decay, but this is not necessarily true for
the states with binding energies $9/n^{2}$. A state with $L^{\prime}=L-1$ will
be available if $9/n^{2}<6.25/m^{2}$ for some $m$ with $L\leq m<n$, i.e. if
$m<n\sqrt{6.25/9}=n/1.2$ is consistent with $m\geq L$. This will be true for
$n>1.2L$, so the dangerous range is states with $L<n\leq 1.2L$. In order for
this range to include an integer, we must have $L>5$. For example, consider
the spin-triplet state with $L=7$ and $n=8$, with dimensionless binding energy
$9/8^{2}\simeq 0.14$ in the Coulombic limit. The lowest-lying $L=6$ spin-
triplet state that is accessible via $\Delta L=1$, $\Delta S=0$ transitions
has $n=7$, and consequently binding energy $6.25/7^{2}=0.13$ in the Coulombic
limit; thus the $L=7$ state cannot decay through such a transition. For the
$L=5$ case, in the Coulombic limit the states are degenerate, and so we would
need to perform a more careful calculation.

If we now consider the case of even-$L$ states, the situation is reversed
between the spin-singlet and the spin-triplet; in the spin-triplet case we
expect the even-$L$ states will always be able to decay to their accompanying,
more deeply bound state with $L^{\prime}=L-1$, with the exception of the $L=0$
case where no such state exists (we will consider the $L=0$ case below). In
the spin-singlet case, the same argument as previously tells us that for
$L>0$, a state with $L^{\prime}=L-1$ will be available except in the case
where $L<n\leq 1.2L$, which in the case of even $L$ is potentially relevant
for $L\geq 6$.

Therefore, based on the Coulombic limit, we would predict that the only
possible (meta)stable states with $L>0$ are $L+S$-even states (spin-triplet
with $L$ odd or spin-singlet with $L$ even) with $L\geq 5$, $L<n\leq 1.2L$,
and the eigenstate structure corresponding to the $Z=6$ eigenvalue. These
high-$L$ states may not even be bound for masses of interest to us, and in any
case the capture rate into them is likely to be very small.

###  F.2 Results for general representations

Let us now consider the more general situation where the DM is the lightest
component of an SU(2) multiplet in a real representation of odd dimension $N$.
Larger representations require higher DM masses to obtain the correct relic
density (e.g. Ref. Bottaro:2021snn ), and hence the unbroken-SU(2)
approximation is likely to be better at their thermal masses. Recall, however,
that the condition to be in the adiabatic regime is mass-independent,
$v\lesssim\delta/m_{\scriptscriptstyle W}$, so while this is a feature of the
broken SU(2) symmetry, we expect it to be retained at sufficiently low
velocities even for very heavy DM.

In the unbroken-SU(2) limit we can use the results of Ref. Mitridate:2017izz
for general representations. They proceed by decomposing the two-particle
state into eigenstates of isospin $I$ ($I=1,3,\cdots,2N-1$); this corresponds
to identifying the eigenstates of the potential in our language. They find the
eigenvalue associated with the state with isospin $I$ is
$\lambda=(2N^{2}-1-I^{2})/8$ (where positive eigenvalues correspond to
attracted states, as per our previous convention) and so the most-attracted
channel is the singlet, where $\lambda=(N^{2}-1)/4$. The isospin singlet
corresponds to an $L+S$-even state with total charge $Q=0$, and for the
quintuplet $\lambda=6$, as discussed above. In general, states with
$I<\sqrt{2N^{2}-1}$ can support bound states; for the quintuplet this means we
have $I=1,3,5$ bound states.

In the adiabatic regime, SU(2)-breaking effects cause the lowest-energy state
at large distances (i.e. the $L+S$-even, $Q=0$ state of two identical DM
particles) to smoothly transition into the isospin-singlet state at short
distances. Thus in this regime, we expect both the Sommerfeld-enhanced direct
annihilation and the bound state capture rates to be consistent with an
initial isospin-singlet state. Since the dominant bound-state capture process
(dipole emission of a gauge boson) changes isospin by 2, the final state must
then have $I=3$, i.e. it is a SU(2) adjoint (this state is $L+S$-odd and the
three components have $Q=0,\pm 1$).

Consequently, within the adiabatic approximation, we only need to concern
ourselves with (Sommerfeld enhanced) direct annihilation from the isospin-
singlet state, and capture from an isospin-singlet initial state to an
isospin-triplet final state, with the relevant eigenvalues being
$\lambda_{i}=(N^{2}-1)/4$ and $\lambda_{f}=(N^{2}-5)/4$. (Transitions amongst
bound states can involve higher-$I$ states; in particular $I=5$ for the
quintuplet, with $\lambda=3$.)

We will thus focus in this appendix on singlet-to-adjoint transitions. We
reiterate that this approximation is not appropriate if the gauge symmetry is
actually unbroken or at the high velocities associated with thermal freezeout
in the early universe (as studied in e.g. Ref. Binder:2023ckj ), where other
transitions can also contribute significantly and may dominate. The quality of
this approximation—i.e. the degree to which the incoming state retains non-
singlet components at small $r$, which could contribute significantly to the
capture rate—is an interesting question, but we ignore it here, as our main
purpose is simply to develop some simple intuition for the importance of
bound-state effects for the gamma-ray endpoint signal. The corresponding
approximation for the quintuplet appears to do a reasonable job of estimating
the relative size of bound-state capture and annihilation, as we see in Fig.
4.

For these singlet-to-adjoint transitions, we can write the group theory
coefficients for bound state formation from Ref. Mitridate:2017izz  in the
simplified form:

$\displaystyle C^{a1b}_{\mathcal{J}}$
$\displaystyle=\frac{1}{\sqrt{T_{R}d_{R}}}\text{Tr}(T^{b}T^{a}),$ (153)
$\displaystyle C^{a1b}_{\mathcal{\tau}}$
$\displaystyle=i\frac{1}{\sqrt{T_{R}d_{R}}}\text{Tr}(T^{b}T^{c}T^{d})f^{acd}=-\frac{1}{\sqrt{T_{R}d_{R}}}\text{Tr}(T^{b}T^{c}T^{d})(T^{a}_{\text{adj}})^{cd}.$

We can now use $\text{tr}(T^{a}T^{b})=T_{R}\delta^{ab}$, and also

$\displaystyle\text{Tr}(T^{b}T^{c}T^{d})(T^{a}_{\text{adj}})^{cd}$
$\displaystyle=\frac{1}{2}T_{R}T_{\text{adj}}\delta^{ab}.$ (154)

Thus, finally we obtain:

$\displaystyle C^{a1b}_{\mathcal{J}}$
$\displaystyle=\sqrt{\frac{T_{R}}{d_{R}}}\delta^{ab},\quad
C^{a1b}_{\mathcal{\tau}}=-\frac{1}{2}\sqrt{\frac{T_{R}}{d_{R}}}T_{\text{adj}}\delta^{ab},$
(155)

where we will show how the $C^{a1b}_{\mathcal{J}}$ and
$C^{a1b}_{\mathcal{\tau}}$ coefficients enter the bound state capture rate in
Secs. F.2.2 and F.2.3.

Now if $R$ is the representation of size $N$, for SU(2) we have
$T_{R}=N(N^{2}-1)/12$, and so $T_{\text{adj}}=2$, while $d_{R}=N$ (and in
particular $d_{\text{adj}}=3$). Thus for SU(2) we obtain the coefficients:

$\displaystyle C^{a1b}_{\mathcal{J}}$
$\displaystyle=\sqrt{\frac{N^{2}-1}{12}}\delta^{ab},\quad
C^{a1b}_{\mathcal{\tau}}=-C^{a1b}_{\mathcal{J}}.$ (156)

Let us also note that we can now extend the argument given in App. F.1.2 to
general representations. A bound state of isospin $I$ will generally have open
decay channels to states with lower isospin (by 2 units) which are hence more
deeply bound due to the larger eigenvalue $\lambda$. The exception is $I=1$
states, which must decay to $I=3$ states which are more shallowly bound for
the same principal quantum number. For this reason, excited $I=1$ states can
be metastable if they have sufficiently large $L$ that all the $I=3$ states
differing by only one unit in $L$ are more shallowly bound. This can occur for
a general representation if
$L<n\leq\frac{N^{2}-1}{N^{2}-5}L=\left(1+\frac{4}{N^{2}-5}\right)L$; this
range will contain an integer if $L>(N^{2}-5)/4$. Thus the threshold $L$ at
which this effect can occur increases as the representation size goes up.

####  F.2.1 Direct annihilation

In this case, if we can evaluate the tree-level cross section for annihilation
from an isospin-singlet initial state to any desired SM final state, we can
account for the Sommerfeld enhancement by simply multiplying the tree-level
cross section by $S=2\pi\alpha_{\scriptscriptstyle W}\lambda_{i}/v$, in the
low-velocity limit. This cross section is given for Majorana fermion DM by
Ref. Mitridate:2017izz  as:

$\displaystyle(\sigma v)_{\text{tree,I=1}}$
$\displaystyle=\frac{\pi\alpha_{\scriptscriptstyle
W}^{2}}{M_{\chi}^{2}}\frac{T_{R}^{2}d_{\text{adj}}}{d_{R}}{}$
$\displaystyle=\frac{\pi\alpha_{\scriptscriptstyle W}^{2}}{2^{4}\times 3\times
M_{\chi}^{2}}N(N^{2}-1)^{2}.$ (157)

Multiplying by the Sommerfeld factor gives:

$\displaystyle(\sigma v)_{\text{I=1}}$
$\displaystyle=\frac{\pi^{2}\alpha_{\scriptscriptstyle W}^{3}}{2^{5}\times
3\times
M_{\chi}^{2}v}N(N^{2}-1)^{3}\rightarrow\frac{\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}\frac{N^{7}}{96},$ (158)

where in the final step we have assumed $N\gg 1$.

Checking, for the wino and quintuplet this yields:

$\displaystyle(\sigma v)_{\text{I=1}}$
$\displaystyle=\frac{\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}\begin{cases}16,&N=3,\\\ 720,&N=5.\end{cases}$ (159)

For the quintuplet this agrees with our calculation above. This also agrees
with the wino result from Ref. Asadi:2016ybp , accounting for our assumption
that the adiabatic condition holds.

####  F.2.2 Capture to the ground state

At small velocities, from Ref. Mitridate:2017izz  we can read off the low-
velocity bound state capture cross section to the $n=1$ state as:

$\displaystyle(\sigma
v)^{n=1,l=0}_{\text{bsf}}=\frac{\pi\alpha_{\scriptscriptstyle
W}^{2}}{M_{\chi}^{2}}\frac{2S+1}{g_{\chi}^{2}}\frac{2^{11}\pi}{3}\sum_{ab}|C^{a1b}_{\mathcal{J}}+(1/\lambda_{f})C^{a1b}_{\mathcal{\tau}}|^{2}\frac{\lambda_{i}^{3}\alpha_{\scriptscriptstyle
W}}{\lambda_{f}v}e^{-4\lambda_{i}/\lambda_{f}}.$ (160)

This expression involves an average over initial states, with degrees of
freedom denoted by $g_{\chi}$; since we are interested in the case where
$100\%$ of the DM captures from the singlet state, we only need to average
over spin degrees of freedom, so $g_{\chi}=2$. Our initial state must be
$L+S$-even and thus for capture to an $L=0$ final state, it must have $S=1$.
Thus we obtain:

$\displaystyle(\sigma v)^{n=1,l=0}_{\text{bsf}}$
$\displaystyle=\frac{\pi\alpha_{\scriptscriptstyle
W}^{2}}{M_{\chi}^{2}}\frac{3}{4}\frac{2^{11}\pi}{3}d_{\text{adj}}\left(\frac{N^{2}-1}{12}\right)|1-(1/\lambda_{f})|^{2}\frac{\lambda_{i}^{3}\alpha_{\scriptscriptstyle
W}}{\lambda_{f}v}e^{-4\lambda_{i}/\lambda_{f}}{}$
$\displaystyle=\frac{8\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}\frac{(N^{2}-9)^{2}(N^{2}-1)^{4}}{(N^{2}-5)^{3}}e^{-4(N^{2}-1)/(N^{2}-5)}.$
(161)

where $d_{\text{adj}}=3$ counts the number of generators and arises from
$\sum_{ab}|\delta^{ab}|^{2}=d_{\text{adj}}$.

In the limit of $N\gg 3$, which is helpful for comparison against direct
annihilation, we obtain the simplified result:

$\displaystyle(\sigma v)^{n=1,l=0}_{\text{bsf}}$
$\displaystyle\rightarrow\frac{\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}8N^{6}e^{-4}\simeq\frac{\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}\frac{N^{6}}{6.8},\quad N\gg 3.$ (162)

(Note in the quintuplet case this approximate value is larger than the truth
by about a factor of 3; it will be a better approximation for larger $N$.)

####  F.2.3 Capture to the $n=2$ states

Now let us consider capture to the $n=2$, $L=1$ states, as their subsequent
decays and annihilations can give rise to endpoint photons, unlike capture
directly to the $n=1$ state. In this case the final state must have $S=0$ (so
the initial state has $L+S$ even). From Ref. Mitridate:2017izz  in the low-
velocity limit the $s$-wave (1st line) and $d$-wave (2nd line) contributions
are:

$\displaystyle(\sigma v)^{n=2,l=1}_{\text{bsf}}$
$\displaystyle=\frac{\pi\alpha_{\scriptscriptstyle
W}^{2}}{M_{\chi}^{2}}\left(\frac{2S+1}{g_{\chi}^{2}}\right)\frac{2^{12}\pi\lambda_{i}}{9\lambda_{f}^{5}}\frac{\alpha_{\scriptscriptstyle
W}}{v}\left[\sum_{ab}|C^{a1b}_{\mathcal{J}}(\lambda_{f}\lambda_{i}(3\lambda_{f}-4\lambda_{i})+C^{a1b}_{\mathcal{\tau}}(-3\lambda_{f}^{2}+12\lambda_{f}\lambda_{i}-8\lambda_{i}^{2})|^{2}\right.{}$
$\displaystyle\left.+\sum_{ab}2^{5}\lambda_{i}^{4}|C^{a1b}_{\mathcal{J}}\lambda_{f}+2C^{a1b}_{\mathcal{\tau}}|^{2}\right]e^{-8\lambda_{i}/\lambda_{f}}{}$
$\displaystyle=\frac{\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}\frac{2^{8}\lambda_{i}}{9\lambda_{f}^{5}}(N^{2}-1)\left[|(\lambda_{f}\lambda_{i}(3\lambda_{f}-4\lambda_{i})-(-3\lambda_{f}^{2}+12\lambda_{f}\lambda_{i}-8\lambda_{i}^{2})|^{2}\right.{}$
$\displaystyle\left.+2^{5}\lambda_{i}^{4}|\lambda_{f}-2|^{2}\right]e^{-8\lambda_{i}/\lambda_{f}}{}$
$\displaystyle=\frac{\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}\frac{2^{4}(N^{2}-1)^{2}}{9(N^{2}-5)^{5}}\left[(N^{6}+9N^{4}-165N^{2}-37)^{2}\right.{}$
$\displaystyle\left.+2^{5}(N^{2}-1)^{4}(N^{2}-13)^{2}\right]e^{-8(N^{2}-1)/(N^{2}-5)}.$
(163)

The first line here agrees with the $s\rightarrow p$ quintuplet result
computed in Eq. (148), once we multiply that result (which was for photon-
mediated capture into a specific $n=2$ $L=1$ bound state) by a factor of 3 to
account for $W$\- and $Z$-mediated capture and a second factor of 3 to account
for the $m=0,\pm 1$ states.

In the limit of large $N$, this expression has the scaling:

$\displaystyle(\sigma v)^{n=2,l=1}_{\text{bsf}}$
$\displaystyle\rightarrow\frac{\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}\frac{2^{4}N^{6}}{9}\left[1+2^{5}\right]e^{-8}\simeq\frac{\pi^{2}\alpha_{\scriptscriptstyle
W}^{3}}{M_{\chi}^{2}v}\frac{N^{6}}{1700}\left[1+32\right],\quad N\gg 5$ (164)

So we see that compared to direct annihilation, in the large-$N$ limit we
expect the various contributions to scale as:

  * •

$p\rightarrow s$ capture to the $n=1,L=0,S=1$ state (contributions to endpoint
photons are power-suppressed): direct annihilation rate $\times 14/N$, in
addition to (when SU(2) is broken) any kinematic suppression of $W/Z$ emission
(by a factor as small as $\alpha/(3\alpha_{\scriptscriptstyle W})$) or
velocity suppression due to the $p$-wave initial state (parametrically
$\mathcal{O}(M_{\chi}^{2}v^{2}/m_{\scriptscriptstyle W}^{2})$ for $v\lesssim
m_{\scriptscriptstyle W}/M_{\chi}$),282828As mentioned above, the large-$N$
approximation overestimates this ratio for the quintuplet by about a factor of
3, and so the rates are actually very comparable.

  * •

$s\rightarrow p$ capture to the $n=2,L=1$ states collectively: direct
annihilation rate $\times 0.06/N$, in addition to any kinematic suppression of
$W/Z$ emission (up to a factor of $3\alpha_{\scriptscriptstyle W}/\alpha$),

  * •

$d\rightarrow p$ capture to the $n=2,L=1$ states collectively: direct
annihilation rate $\times 1.8/N$, in addition to any kinematic suppression of
$W/Z$ emission (up to a factor of $3\alpha_{\scriptscriptstyle W}/\alpha$) or
velocity suppression due to the $d$-wave initial state (parametrically
$\mathcal{O}(M_{\chi}^{4}v^{4}/m_{\scriptscriptstyle W}^{4})$ for $v\lesssim
m_{\scriptscriptstyle W}/M_{\chi}$).

We see that the only contribution that is not suppressed at low velocities and
that gives rise to leading-power contributions to the endpoint photon spectra
(via its decay and subsequent annihilation) is generically expected to have a
cross section 2 or more orders of magnitude below direct annihilation.
Consequently, it is quite plausible for bound state formation to be a large or
even dominant contribution to the inclusive annihilation rate when the
velocity suppression for higher partial waves is mild or absent and the gauge
bosons are massless (as in the case of freezeout), while simultaneously having
a generically small effect on the endpoint spectrum for indirect detection,
especially at low velocities ($v\ll m_{\scriptscriptstyle W}/M_{\chi}$) or
where the $n=2$ states are too loosely bound to allow $W$\- or $Z$-mediated
capture.

One might ask about the contribution from capture into states with $n>2$.
While in the unbroken limit bound states with large $L$ and $n$ may play a
large role in the capture rate (e.g. Binder:2023ckj ), for the parameter space
we have considered in this paper, the number of bound states is always
truncated by the non-zero SU(2) breaking scale, preventing large enhancements
from the proliferation of high-$n$ states. Furthermore, we expect velocity
suppressions (of order $(M_{\chi}v/m_{\scriptscriptstyle W})^{2L}$) for all
capture rates with initial $L>0$. Finally, within our adiabatic approximation
both initial and final states always experience attractive interactions, with
couplings that obey $\lambda_{i}/\lambda_{f}>1$, leading to increasingly
strong exponential suppression for large-$n$ states and avoiding the
potentially unitarity-violating region of parameter space identified in Ref.
Binder:2023ckj .

## References

  * (1) M. Cirelli, N. Fornengo, and A. Strumia, Minimal dark matter, Nucl. Phys. B753 (2006) 178–194, [hep-ph/0512090].
  * (2) A. Mitridate, M. Redi, J. Smirnov, and A. Strumia, Cosmological Implications of Dark Matter Bound States, JCAP 05 (2017) 006, [arXiv:1702.01141].
  * (3) S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo, and L. Vittorio, Closing the window on WIMP Dark Matter, Eur. Phys. J. C 82 (2022), no. 1 31, [arXiv:2107.09688].
  * (4) M. Cirelli, A. Strumia, and M. Tamburini, Cosmology and Astrophysics of Minimal Dark Matter, Nucl. Phys. B787 (2007) 152–175, [arXiv:0706.4071].
  * (5) M. Beneke, R. Szafron, and K. Urban, Wino potential and Sommerfeld effect at NLO, Phys. Lett. B 800 (2020) 135112, [arXiv:1909.04584].
  * (6) K. Urban, NLO electroweak potentials for minimal dark matter and beyond, JHEP 10 (2021) 136, [arXiv:2108.07285].
  * (7) F. Aharonian, D. Khangulyan, and D. Malyshev, Cold ultrarelativistic pulsar winds as potential sources of galactic gamma-ray lines above 100 GeV, Astron. Astrophys. 547 (2012) A114, [arXiv:1207.0458].
  * (8) G. Sinnis, A. Smith, and J. E. McEnery, HAWC: A Next generation all - sky VHE gamma-ray telescope, in On recent developments in theoretical and experimental general relativity, gravitation, and relativistic field theories. Proceedings, 10th Marcel Grossmann Meeting, MG10, Rio de Janeiro, Brazil, July 20-26, 2003. Pt. A-C, pp. 1068–1088, 2004.  astro-ph/0403096.
  * (9) HAWC Collaboration, A. U. Abeysekara et al., Sensitivity of HAWC to high-mass dark matter annihilations, Phys. Rev. D 90 (2014), no. 12 122002, [arXiv:1405.1730].
  * (10) HAWC Collaboration, A. Albert et al., Dark Matter Limits From Dwarf Spheroidal Galaxies with The HAWC Gamma-Ray Observatory, Astrophys. J. 853 (2018), no. 2 154, [arXiv:1706.01277].
  * (11) HAWC Collaboration, A. U. Abeysekara et al., A Search for Dark Matter in the Galactic Halo with HAWC, JCAP 02 (2018) 049, [arXiv:1710.10288].
  * (12) HAWC Collaboration, A. Albert et al., Search for Dark Matter Gamma-ray Emission from the Andromeda Galaxy with the High-Altitude Water Cherenkov Observatory, JCAP 06 (2018) 043, [arXiv:1804.00628]. [Erratum: JCAP 04, E01 (2019)].
  * (13) HAWC Collaboration, A. Albert et al., Search for gamma-ray spectral lines from dark matter annihilation in dwarf galaxies with the High-Altitude Water Cherenkov observatory, Phys. Rev. D 101 (2020), no. 10 103001, [arXiv:1912.05632].
  * (14) H.E.S.S. Collaboration, J. A. Hinton, The Status of the H.E.S.S. project, New Astron. Rev. 48 (2004) 331–337, [astro-ph/0403052].
  * (15) H.E.S.S. Collaboration, F. Aharonian et al., H.E.S.S. observations of the Galactic Center region and their possible dark matter interpretation, Phys. Rev. Lett. 97 (2006) 221102, [astro-ph/0610509]. [Erratum: Phys.Rev.Lett. 97, 249901 (2006)].
  * (16) H.E.S.S. Collaboration, F. Aharonian et al., Observations of the Sagittarius Dwarf galaxy by the H.E.S.S. experiment and search for a Dark Matter signal, Astropart. Phys. 29 (2008) 55–62, [arXiv:0711.2369]. [Erratum: Astropart.Phys. 33, 274–275 (2010)].
  * (17) H.E.S.S. Collaboration, A. Abramowski et al., Search for a Dark Matter annihilation signal from the Galactic Center halo with H.E.S.S, Phys. Rev. Lett. 106 (2011) 161301, [arXiv:1103.3266].
  * (18) H.E.S.S. Collaboration, A. Abramowski et al., Search for Photon-Linelike Signatures from Dark Matter Annihilations with H.E.S.S., Phys. Rev. Lett. 110 (2013) 041301, [arXiv:1301.1173].
  * (19) H.E.S.S. Collaboration, A. Abramowski et al., Search for dark matter annihilation signatures in H.E.S.S. observations of Dwarf Spheroidal Galaxies, Phys. Rev. D 90 (2014) 112012, [arXiv:1410.2589].
  * (20) H.E.S.S. Collaboration, A. Abramowski et al., Constraints on an Annihilation Signal from a Core of Constant Dark Matter Density around the Milky Way Center with H.E.S.S., Phys. Rev. Lett. 114 (2015), no. 8 081301, [arXiv:1502.03244].
  * (21) H.E.S.S. Collaboration, H. Abdallah et al., Search for dark matter annihilations towards the inner Galactic halo from 10 years of observations with H.E.S.S, Phys. Rev. Lett. 117 (2016), no. 11 111301, [arXiv:1607.08142].
  * (22) H.E.S.S. Collaboration, H. Abdalla et al., H.E.S.S. Limits on Linelike Dark Matter Signatures in the 100 GeV to 2 TeV Energy Range Close to the Galactic Center, Phys. Rev. Lett. 117 (2016), no. 15 151302, [arXiv:1609.08091].
  * (23) H.E.S.S. Collaboration, H. Abdallah et al., Search for $\gamma$-ray line signals from dark matter annihilations in the inner Galactic halo from ten years of observations with H.E.S.S, Phys. Rev. Lett. 120 (2018) 201101, [arXiv:1805.05741].
  * (24) H.E.S.S. Collaboration, H. Abdalla et al., Searches for gamma-ray lines and ’pure WIMP’ spectra from Dark Matter annihilations in dwarf galaxies with H.E.S.S, JCAP 11 (2018) 037, [arXiv:1810.00995].
  * (25) H.E.S.S. Collaboration, H. Abdalla et al., Search for Dark Matter Annihilation Signals in the H.E.S.S. Inner Galaxy Survey, Phys. Rev. Lett. 129 (2022), no. 11 111101, [arXiv:2207.10471].
  * (26) J. Holder et al., Status of the VERITAS Observatory, AIP Conf. Proc. 1085 (2009), no. 1 657–660, [arXiv:0810.0474].
  * (27) VERITAS Collaboration, T. Arlen et al., Constraints on Cosmic Rays, Magnetic Fields, and Dark Matter from Gamma-Ray Observations of the Coma Cluster of Galaxies with VERITAS and Fermi, Astrophys. J. 757 (2012) 123, [arXiv:1208.0676].
  * (28) VERITAS Collaboration, S. Archambault et al., Dark Matter Constraints from a Joint Analysis of Dwarf Spheroidal Galaxy Observations with VERITAS, Phys. Rev. D 95 (2017), no. 8 082001, [arXiv:1703.04937].
  * (29) A. Acharyya et al., Search for Ultraheavy Dark Matter from Observations of Dwarf Spheroidal Galaxies with VERITAS, Astrophys. J. 945 (2023), no. 2 101, [arXiv:2302.08784].
  * (30) MAGIC Collaboration, E. Lorenz, Status of the 17-m MAGIC telescope, New Astron. Rev. 48 (2004) 339–344.
  * (31) MAGIC Collaboration, J. Aleksic et al., Searches for Dark Matter annihilation signatures in the Segue 1 satellite galaxy with the MAGIC-I telescope, JCAP 06 (2011) 035, [arXiv:1103.0477].
  * (32) MAGIC, Fermi-LAT Collaboration, M. L. Ahnen et al., Limits to Dark Matter Annihilation Cross-Section from a Combined Analysis of MAGIC and Fermi-LAT Observations of Dwarf Satellite Galaxies, JCAP 02 (2016) 039, [arXiv:1601.06590].