The Quintuplet Annihilation Spectrum

We extend the Effective Field Theory of Heavy Dark Matter to arbitrary odd representations of SU(2) and incorporate the effects of bound states. This formalism is then deployed to compute the gamma-ray spectrum for a 5 of SU(2): quintuplet dark matter. Except at isolated values of the quintuplet mass, the bound state contribution to hard photons with energy near the dark-matter mass is at the level of a few percent compared to that from direct annihilation. Further, compared to smaller representations, such as the triplet wino, the quintuplet can exhibit a strong variation in the shape of the spectrum as a function of mass. Using our results, we forecast the fate of the thermal quintuplet, which has a mass of $\sim$13.6 TeV. We find that existing H.E.S.S. data should be able to significantly test the scenario, however, the final word on this canonical model of minimal dark matter will likely be left to the Cherenkov Telescope Array (CTA).


Introduction
A fundamental difficulty in the search for the particle nature of dark matter (DM) is the enormous array of well motivated candidates that exist. There are many ways to forge a path through the vast model space. One guiding principle is a bottom-up notion of simplicity-a preference is given for the minimal modification to the Standard Model (SM) consistent with observations. As originally emphasized in Ref. [1], with the mantra of minimality, few models are simpler than quintuplet DM. Within this model, the SM is augmented by a single new field, a Majorana fermion transforming in the 5 or quintuplet representation of SU (2), and as a gauge singlet under SU(3)×U(1). This representation sits at a "sweet spot." It is large enough that no additional symmetries are needed to make it cosmologically stable, with decay operators only appearing at dimension-6. Simultaneously, it is small enough that the SU(2) Landau pole remains above the GUT scale. After electroweak symmetry breaking, the five Majorana fermions of the quintuplet reorganize themselves into a neutral Majorana fermion, χ 0 -the DM candidate -as well as singly and doubly-charged Dirac fermions, χ + and χ ++ , which will play important roles in the phenomenology. The coupling between the DM and the SM is fixed by measurements of the electroweak coupling α W , leaving the DM mass, M χ , as the single free parameter in the model. Assuming a conventional thermal origin of the DM, even the mass becomes fixed to M χ = 13.6 ± 0.8 TeV [2,3]. 1 However, if the quintuplet represents only a fraction of the DM or was produced through a non-standard cosmology, a wider range of TeV-scale masses becomes possible. A brief review of quintuplet DM and our conventions for it is provided in App. A. The question of whether the quintuplet is the DM of our Universe can be probed through indirect detection searches for its annihilation signal, χ 0 χ 0 → γ + X. Here, γ represents the experimentally detected hard photon with energy comparable to M χ , and X represents additional undetected radiation. In this context, we consider a hard photon to be one with energy near the mass of the DM, so E γ ∼ M χ , and these photons are of particular interest as they give a line-like signature in the photon spectrum at the multi-TeV scale, which is not expected from any plausible astrophysical backgrounds. 2 There are reasons to be optimistic that the coming decade may bring with it the detection of such TeV photons from DM annihilation, primarily on account of the rapid progress in the experimental program of searching for TeV-PeV photons. Building on the successes of DM searches with ongoing air and water Cherenkov telescopes such as HAWC [8][9][10][11][12][13], H.E.S.S. [14][15][16][17][18][19][20][21][22][23][24][25], VERITAS [26][27][28][29], 1 This value includes the contribution of bound states to the relic abundance. The value without these effects is closer to 9.6 TeV [4]. The value is, however, computed with the LO rather than NLO potentials of Refs. [5,6] that we will make use of throughout, as described in Sec. 3.1. 2 As discussed in Ref. [7], astrophysical emission processes will almost exclusively generate spectra that are broader than the energy resolution of instruments such as H.E.S.S. and CTA. The possible exception explored in that work is inverse Compton emission deep in the Klein-Nishina regime. In that case, where the photons acquire essentially the entire electron or positron momentum, the photon spectrum is determined solely by the initial electron and positron distribution. Sufficiently narrow leptonic distributions to generate an apparent gamma-ray line could be produced by pulsar winds, but such a scenario is by no means generic, and the spatial distribution of such a signal would ultimately be distinct to that predicted for DM.
The central focus of the present work is to take this possibility seriously, and therefore derive a precise theoretical prediction for the photon spectrum from quintuplet annihilation. Even though the quintuplet represents at most a single-parameter model characterized solely by its mass, there are a number of effects that must be carefully accounted for in order to produce an accurate spectrum. Arguably the most important of these is the fact that as two χ 0 approach one another, once they come within a distance r ∼ m −1 W , they can exchange virtual electroweak bosons, and thereby experience a potential that can significantly perturb their initial wavefunctions. This effect is referred to as Sommerfeld enhancement, and can modify the expected annihilation rate by orders of magnitude [4,[48][49][50][51]. Further, the problem involves several hierarchies of scale, which necessitates an effective field theory treatment. The large hierarchy between the DM mass and the electroweak scale, M χ ≫ m W , manifests itself through large Sudakov double logarithms, α W ln 2 (M χ /m W ), which lead to a breakdown in naive perturbation theory. Perturbative control can be restored by resumming these logarithms using the techniques of effective field theory (EFT), in particular one built using Soft-Collinear Effective Theory (SCET) [52][53][54][55]. This solution has previously been implemented in heavy dark-matter models of neutralinos, including the wino and higgsino [56][57][58][59][60][61][62][63][64][65][66][67]. An additional source of large logs occurs due to our insistence on hard photons near the endpoint of the spectrum, with E γ ∼ M χ . More carefully, we focus on photons with z = E γ /M χ , such that (1−z) ≪ 1. This kinematic restriction gives rise to large terms of the form ln(1−z), which can also be resummed as shown in, for instance, Ref. [62]. In summary, the effects of Sommerfeld enhancement, electroweak Sudakov logarithms, and endpoint contributions have been incorporated in a number of scenarios of neutralino DM. For instance, Refs. [62,64] developed an EFT framework for including all of these effects in the case of wino-like DM, where DM is a 3 or triplet of SU (2), allowing the spectrum to be computed to next-to-leading logarithmic (NLL) accuracy. However, the power of EFT and factorization is that many aspects of those calculations should not depend on the specific DM representation, so that results for the wino -an SU(2) triplet -should lift to the quintuplet and in fact, as we will discuss, to any higher representation. In Sec. 2 of the present work we will demonstrate this explicitly, and in particular compute the NLL quintuplet direct annihilation spectrum.
The quintuplet also represents an opportunity to extend this formalism to incorporate an additional (potentially) important source of photons: bound state formation and decay. The starting point for these bound states is similar to that of the Sommerfeld enhancement to annihilation. As there, once the initial χ 0 pair comes within r ∼ m −1 W , they can exchange a W -boson and thereby convert into a χ + χ − pair. This pair of charged particles can now emit a photon, capturing into a bound state, which can be thought of as the DM quintuplet analogue of positronium. At higher DM masses, it is also possible to emit a W boson and capture into a charged bound state, and further on-shell Z boson emission becomes a significant contributor to the formation of neutral bound states. All of these bound states are generally unstable and may decay either to lighter bound states in the spectrum, or directly to SM particles. In principle the contribution of bound states could also be important for the wino, however, as shown in Ref. [68] it is irrelevant for present-day indirect detection (albeit possibly relevant in the early Universe, see for example Ref. [2]). For the thermal wino, there is only one bound state present in the spectrum, and capture to it via emission of a dipole photon is forbidden by spin statistics. For heavier winos or wino-like particles (with non-thermal histories to prevent overclosure of the Universe), there is a non-zero rate for capture to a spectrum of bound states via dipole photon emission (and at sufficiently high masses, W and Z emission). However, dimensionless numerical factors arising from the wavefunction overlap integral render these rates small compared to direct annihilation, even in the limit of very heavy DM where SU (2) can be treated as approximately unbroken.
For the quintuplet, it was already anticipated in Ref. [68] that the higher preferred mass would allow for a more complicated spectrum of bound states, and the suppression from dimensionless factors should not be as severe. Furthermore, previous studies of quintuplet DM have found that bound state formation plays a key role in setting the abundance of DM in the early Universe [2] and can be non-negligible in indirect searches for quintuplet annihilation in the Milky Way halo [69]. In the present work, we apply our SCET-based formalism to precisely compute contributions to the annihilation signal from bound state formation followed by decay. The existence of bound states in the quintuplet spectrum roughly induces the following independent (i.e. non-interfering) contribution to the annihilation spectrum, where the sum is over the set of bound states in the theory, σ(χ 0 χ 0 → B + X us ) is the production cross section for the bound state which happens via the emission of an unobserved ultra-soft particle X us = γ, W, Z, Γ B is the total decay width for the bound state, and dΓ B→γ+X /dz is the differential decay rate to a measured photon carrying energy zM χ . This result is essentially a consequence of the narrow width approximation; we will justify it in App. C. Taking Eq. (1.1) as given, the problem of introducing the bound state contribution is reduced to computing three ingredients. The first of these is the capture cross section, σ(χ 0 χ 0 → B + X us ), or heuristically the probability for a given bound state to form. We compute this using first-order perturbation theory in non-relativistic quantum mechanics, following Refs. [68,70], with an extension to handle bound states produced by W/Z emission instead of photon emission. The second of these is determining the fully inclusive decay rate for the bound state, Γ B . To obtain this result we need both the decay rates for excited states to transition to lower-energy states (calculated using perturbation theory, similar to the initial capture process) and the rates for decay via annihilation into SM particles. The final ingredi-ent is to compute the photon spectrum of the bound state decay. As emphasized above we are primarily interested in the spectrum of hard photons, with E γ ∼ M B /2. This implies that the photons are kinematically restricted to the endpoint region, and this will allow us to draw on the SCET-based machinery for calculating endpoint spectra developed in Refs. [62,64]. For this final step, we will also need to consider additional operators mediating the hard process that are not needed for the direct annihilation process. The reason for this is that the decaying bound states may have different spin and angular momentum quantum numbers than those configurations that dominate the direct annihilation. Our calculation of the aforementioned bound state effects is divided between two sections: in Sec. 3 we compute the formation rates, while in Sec. 4 we determine the fate of each bound state, accounting for transitions between bound states and their ultimate decay to SM particles. When these results are combined, we observe that the contribution of the bound state annihilation to the hard photon spectrum is only a few percent of that from direct annihilation at most masses. The contribution from bound state annihilation is small for a number of reasons: firstly, bound-state formation from a χ 0 χ 0 initial state with L = 0 requires capture into an excited state, which is generically suppressed by a wavefunction overlap factor, and the transition from the L = 0 state to the lowest-lying L = 1 states also has an accidentally small numerical prefactor in the cross section-for the wino, this coefficient is zero in the limit of unbroken SU (2). If we instead consider L > 0 initial states, these contributions are velocity suppressed due to the small non-relativistic speed of DM in our Milky Way Galaxy. In addition, odd-L initial states must have S = 1 to ensure the asymmetry of the χ 0 χ 0 wavefunction. They thus give rise to S = 1 bound states (ignoring spin-flip transitions, which are suppressed); we find that such bound states have power-suppressed contributions to the endpoint photon spectrum. All these effects are discussed in App. F, where we also show that in the limit of high DM mass and large representation size, the ratio of bound-state formation to direct annihilation is expected to decrease further for larger representations (in the context of indirect detection of hard gamma rays).
In Sec. 5 we combine the pieces to obtain the full quintuplet endpoint spectrum and annihilation cross section as a function of DM mass. As for the wino, the annihilation cross section exhibits a rich structure and rapid variation associated with near-zero energy bound states, characteristic of the Sommerfeld enhancement. Unlike the wino, however, we also see a strong variation in the shape of the spectrum itself: the energy distributions of photons resulting from a quintuplet annihilation can depend sensitively on its exact mass. Using these results, we estimate the sensitivity of existing H.E.S.S. data to the thermal quintuplet, finding that for commonly adopted DM profiles in the Milky Way, the signal should already be either observable or in tension. If we adopt a more conservative DM profile, however, our estimate is that the final word on the quintuplet will await the data that will soon be collected by CTA. Finally, our conclusions are presented in Sec. 6, with several extended discussions and details relegated to appendices. TeV, the resummed cross-section is more than a factor of four smaller. (Right) The cross section for the thermal quintuplet integrated over z ∈ [z cut , 1]. A significant reduction in the theoretical uncertainty is obtained for the NLL result. Further, we see an increase in the cross section in the range appropriate for the H.E.S.S. energy resolution, which arises due to the presence of endpoint photons. (The right figure can be directly contrasted to that of the wino, presented in Ref. [64].) Both figures consider solely the contribution from direct annihilation; incorporating the bound state contributions will be a major focus of this work.

Direct Annihilation
While bound states represent a novel addition to the spectrum for the case of the quintuplet, direct annihilation via χ 0 χ 0 → γ + X remains the dominant contribution to the hard photon spectrum for most masses, and we will compute it in this section. To do so we will draw on the EFT formalism developed to compute the leading log (LL) spectrum in Ref. [62], and then extended to NLL in Ref. [64]. The formalism there was applied to the wino, yet as emphasized in those references the approach can be readily extended to other DM candidates, especially to cases where the DM is simply charged under SU(2). This section represents an explicit demonstration of that claim. We begin in Sec. 2.1 by briefly reviewing the framework developed in Refs. [62,64]. In doing so, we will focus on the aspects of the formalism that we will generalize to make it clear how to extend the calculation to additional SU(2) representations of DM, and we defer to those references for a complete discussion of all the relevant ingredients.
Having done this, in Sec. 2.2 we will then demonstrate explicitly how the calculation can be extended to the quintuplet, and provide results for the LL and NLL spectrum. Before moving into the details, however, let us already demonstrate the importance of performing the NLL resummation, and the precision obtained by doing so. In the left of Fig. 1 we show a comparison of the line cross section computed to NLL accuracy, with the associated uncertainties, compared to the result if we had only performed a tree level computation of the rates, augmented with the Sommerfeld enhancement determined from the NLO potentials of Ref. [71] (discussed further in Sec. 3.1). Note, the line cross section is the annihilation rate to a two-photon final state, specifically the rate for χ 0 χ 0 → γγ + half the rate for χ 0 χ 0 → γZ, with the photons having an energy E γ = M χ . We see that at larger masses the difference between the two methods can be significant. The NLL result is already a factor of four smaller than the tree level approximation at the thermal mass, and by 100 TeV the difference is more than an order of magnitude. The EFT formalism also incorporates endpoint photons (described in detail below) that have E γ ≲ M χ , which, given the finite energy resolution of IACTs like H.E.S.S., are indistinguishable from the line. This is shown on the right of Fig. 1, where we plot the integrated cumulative cross section down to a given z cut , which is almost a factor of two larger within the H.E.S.S. resolution as compared to the value for z cut = 1. In addition, this figure demonstrates the considerable reduction in theoretical uncertainty of the cross section obtained at NLL.

Review: An EFT for the endpoint spectrum from DM annihilation
As described in the introduction, a precise prediction of the hard photon spectrum arising from the annihilation of heavy DM charged under SU(2) mandates an accounting of a number of physical effects. The benefit of approaching this problem using the EFT formalism reviewed in this section is that the various effects will factorize; heuristically, we will be able to separate the physics associated with different scales into objects that can be calculated independently.
To reiterate, the calculation we are interested in is the spectrum, dσ/dz, of hard photons resulting from DM annihilations, where "hard" means that we are interested in photons carrying away large energy fractions z = E γ /M χ ∼ 1. The starting point is two incoming neutral DM particles, χ 0 , which are asymptotically momentum eigenstates described by plane waves. These states, initially with momenta p 1,2 ∼ M χ v ∼ 10 −3 M χ (in the Milky Way's halo) will eventually be within a distance r ∼ m −1 W , at which point they will experience an interaction potential that perturbs their wave-functions away from plane waves. The perturbed wavefunctions can yield a significantly enhanced probability for the particles to have a separation r ∼ M −1 χ , where the hard annihilation occurs, and thereby provide a large boost to the cross section. Restricting our attention to the case where M χ ≫ m W , the Sommerfeld enhancement occurs on a parametrically larger distance scale than the annihilation, and so we can factorize it out from the cross section as follows, 3 Here a ′ b ′ ab are adjoint SU(2) indices, and for the wino we can write Here χ v is the field describing the non-relativistic DM, and the label v appears just as for a heavy quark in Heavy Quark Effective Theory (HQET). The non-relativistic DM effective theory that governs the dynamics of the field is reviewed in Ref. [62], and as shown there the above expressions for F a ′ b ′ ab can be directly related to conventional Sommerfeld enhancement factors, as we now review. In the broken phase, the triplet wino is described by a neutral Majorana fermion χ 0 , and a heavier charged Dirac fermion χ ± . DM annihilations proceed through the neutral states, and to ensure the antisymmetry of the state, at lowest order in the DM velocity (s-wave), the initial state must be a spin singlet, which we represent through the notation (χ 0 χ 0 ) S . From this initial state, Eq. (2.2) describes the fact that through the exchange of electroweak bosons, not only will the incident wave-functions be perturbed, but further there are two final states that the system could evolve into by the time the hard process is initiated, χ 0 χ 0 or χ + χ − . In Ref. [62] these matrix elements were determined as, Here s 00 and s 0± are the Sommerfeld factors that need to be computed, and note that if the Sommerfeld effect were neglected, they would take the values s 00 = 1 and s 0± = 0. We emphasize again that the above expressions only hold for the wino; for the quintuplet, we would also need to account for the presence of the doubly-charged states.
To NLL accuracy, with the Sommerfeld enhancement stripped at the first stage of the matching, the differential cross section with factorized dynamics in SCET is given in terms of the hard function H, the jet functions J ′ n , J γ , and the soft function S ′ [62,64] dσ dz The cross section can then be refactorized into a combination of the following seven factors [62,64] dσ dz where we have suppressed the color indices on the hard scattering cross section in Eq. (2.1). Let us briefly provide some intuition for this expression. The DM annihilation will occur an astrophysical distance from our telescope, and therefore no matter how complex the final state is, we should only expect to see a single photon from the decay. This implies we are sensitive to the annihilation χ 0 χ 0 → γ + X, where we must be inclusive over the unobserved X. Although unobserved, X cannot in fact be completely arbitrary. Our choice to search for photons with energy E γ = zM χ , with (1 − z) ≪ 1, implies that the invariant mass of the recoiling states is constrained to be small: m X = 2M χ √ 1 − z ≪ M χ . This implies that the spray of radiation the photon recoils against must be a jet. With this picture in mind, we can apply a conventional SCET factorization to our problem, breaking it into a function for the hard scattering (H), the collinear radiation in the direction of the photon (J γ ) and the recoiling jet (J ′ n ), and finally the soft wide angle radiation (S ′ ). As shown in Ref. [62], this factorization can be achieved, but it is insufficient to fully separate the scales that appear when accounting for the finite masses of the electroweak bosons. For this, one must further factorize J ′ n into the two functions Jn and H Jn , and S ′ into the three functions S, H S , and C S . The full details of this argument, together with the field theoretic definition of each function, is provided in Ref. [62], with Ref. [64] demonstrating that the factorization remains valid even when computing to NLL accuracy. The central utility to Eq. (2.5) is that each of the functions can be computed separately -and independently of the DM representationand then brought to a common scale using renormalization group evolution. This facilitates a full result which resums logarithms of m W /M χ , but also of (1 − z), which can be large given that we are searching for photons near the endpoint with z ∼ 1.
The above factorization is perfectly sufficient for the wino. However, the fact that the DM matrix elements in Eq. (2.2) index the DM with an adjoint SU(2) label demonstrates that this form can only be appropriate for the wino, and further implies that dσ/dz is not representation independent. We will now generalize this. To do so, we must revisit the matching of the full theory to our EFT, as this is where the DM representation enters the calculation. Doing so, it is straightforward to show that the tree-level matching can be achieved by way of a single hard scattering operator, given by with Wilson coefficient, In the above result T χ are the generators of SU (2), and therefore c, d correspond to adjoint indices. However, T χ is written in whatever representation is appropriate for DM; a review of the relevant form for the generators in the triplet and quintuplet representations is given in App. A. Equation (2.6) provides the hard operator before a BPS field redefinition [55]. This transformation must be performed in order to factorize the interactions of the heavy DM from the ultrasoft radiation. Accordingly, we now perform a field redefinition, where S v and Y aa ′ n are both ultrasoft Wilson lines, but the former is in the v direction and the same representation as DM, whereas the latter is in the n direction and adjoint representation.

The operator then transforms as
where in the final line we have defined, 4 To arrive at this result, we made use of the identity To be explicit, the identity has allowed us to replace a pair of S v Wilson lines, which are in the DM representation, with a single adjoint Y v Wilson line.
For the wino, the anticommutator can be readily evaluated, Substituting this into Eq. (2.9), and using with Wilson coefficients exactly matching those determined in Ref. [58]. Nevertheless, in order to fully factorize the DM representation, we should keep the anticommutator unexpanded. In particular, by so doing only the combination χ T v iσ 2 T a χ , T b χ χ v retains any knowledge of the DM representation. We can then introduce a modified definition of Eqs. (2.1) and (2.2) which achieves a complete separation of the DM representation from the factorized expressions in Eq. (2.5). In detail, we write (2.14) With this factorization we move the DM dependence into F χ , and the remaining factors in Eq. (2.9) determine the hard matching onto the SCET calculation of dσ/dz. Importantly, 4 We note that in order to reproduce the operator definitions in Ref. [58] and the works that followed it, , where all indices have been transposed. We believe the index ordering in that work simply had a typo, and note that to the order all wino calculations have been performed so far, flipping these indices would not impact the results.
when the DM factors are stripped, what remains in O is Y abcd B ic ⊥n B jd ⊥n iϵ ijk (n −n) k , which exactly matches the DM stripped contribution in O 2 . This implies that if we intend to compute the cross section for the wino using Eq. (2.14), we need to make two minor changes to the approach used to compute it with Eqs. (2.1) and (2.2). Firstly, we must determine the new DM matrix elements F a ′ b ′ ab χ in terms of the Sommerfeld factors s 00 and s 0± , specified in Eq. (2.3), and then evaluate the contraction of F χ into dσ/dz. Secondly, in the SCET calculation, we match onto the hard function with C 1 (µ) = 0 and C 2 (µ) = −πα W (µ)/2M χ as opposed to the values in Eq. (2.13). Of course, there was nothing fortuitous in the fact that the DM representation can be factored out so simply, this is simply a manifestation of the power of the EFT approach developed in Refs. [62,64]. The factorization in Eq. (2.5) is determined by the relevant degrees of freedom in the theory at scales below M χ , and this is not altered by changing the DM representation, so the factorization remains unaffected.
To make this point explicit, let us demonstrate that at LL these two approaches yield the same result for the wino; a straightforward generalization of the argument below confirms this conclusion persists at NLL. As outlined above, the calculation in Ref. [62] is modified in two ways: a new set of DM matrix elements is computed in F χ , and then in the SCET calculation an alternative matching is provided onto the hard function, H. After the BPS field redefinition, the DM representation contracts into the ultrasoft Wilson lines, and therefore in the SCET calculation is contracted into the soft function,S. (Here,S = C S S is the combination of the collinear-soft and soft functions from Eq. (2.5).) In full, what we will need to recompute is the combination HH SS , as the first and last of these factors is modified, and although H S remains unchanged, it connects these two objects.
Let us begin by reviewing the relevant part of the calculation as it appeared in Ref. [62]. Accounting for the running of the hard function between 2M χ and m W , we have where the renormalization group evolution is encoded in U H , which is defined in terms of α W = α W /4π and C A = 2 the SU(2) adjoint Casimir. The tree level matching coefficients H tree i , are determined by the operator Wilson coefficients C 1 and C 2 , in particular H tree Given the two different structures of the ultrasoft Wilson lines in Eq. (2.12), there are four soft functions at the amplitude square level that evaluate to,S (2.16) These must then be contracted into F a ′ b ′ ab as defined in Eq. (2.2), which using Eq. (2.3) can be evaluated in terms of s 00 and s 0± , with explicit expressions provided in Ref. [62]. The renormalization group evolution of the soft function, and the contraction between it and the hard function, are controlled by H S , which is given by, 5 where U H S quantifies the evolution in a similar fashion to U H . Combining these results, we conclude Re(s 00 s * 0± ) Re(s 00 s * 0± ) (2.18) The functions F 0 and F 1 encode the two combinations of Sommerfeld factors that appeared in parentheses, and appear in the wino LL result. The above is a direct repetition of the calculation performed in Ref. [62], we now demonstrate that the same result is achieved in our modified approach. Firstly, in this approach, the hard matching coefficients are modified, with only H tree 2 non-zero now, because C 1 = 0 (recall, we have a single operator here with the same structure as O 2 ). The only other modifications are the contractions ofS in Eq. (2.16) into F χ rather than F . These can be determined straightforwardly, for instance, The remaining combinations are given by, (2.20) Using these modified results we find that H i H S,ijS in the new basis exactly matches 5 We note that to obtain this result we used H tree S,22 = 1, as opposed to H tree S,22 = 2, which was stated in Ref. [62] that we believe was a typo.
Eq. (2.18), as it must. The utility of this approach, is that having formulated the calculation in this way, if we changed the DM representation, the only part of the calculation that would need to be modified is that the appropriate generalizations of contractions in Eqs. (2.19) and (2.20) would need to be computed. We will evince this by showing that results for the quintuplet can be derived straightforwardly in the next subsection. Again, at NLL an almost identical modification to the approach in Ref. [64] is required, one must simply account for the more complicated forms the hard and soft functions take at that order.

Extension to the quintuplet
In the previous subsection, we reorganized the formalism of Refs. [62,64] in such a way that the dependence on the DM representation is fully encoded in F χ as defined in Eq. (2.14), and explicitly demonstrated this alternative procedure produces the same result at LL. This reorganization has the benefit that the quintuplet calculation (and that for any higher odd-n representation, see e.g. Ref. [3]) is almost identical to that of the wino; there are unique Sommerfeld expressions to compute, and a modification for how the new F χ contracts into the soft Wilson lines in Y abcd , but in essence the computation is the same.
For the Sommerfeld factors, in the broken phase, the five degrees of freedom of the quintuplet reorganize themselves into a neutral Majorana fermion χ 0 , a heavier singly charged Dirac fermion χ ± , and a doubly-charged Dirac fermion χ ±± that is even heavier. (Again, a more complete discussion is provided in App. A.) This implies that there are now three two-body states that can initiate the hard annihilation that are coupled to the initial state through the potential, 6 and we parameterize the various matrix elements as 7 (2.21) Note that if we performed our entire calculation at tree level or neglected the Sommerfeld effect, we would take s 00 = 1 and s 0± = s 0±± = 0. Using these, we can compute the full F χ for an arbitrary set of indices. From these functions we can immediately derive the relevant spectra. At LL, all that is required is to derive the analogue of F 0 and F 1 as they appeared in Eq. (2.18). The main calculation is to compute the equivalent contractions for Eqs. (2.19) and (2.20), which are given bỹ This ignores for the moment the possibility of annihilation through bound states with different quantum numbers, which we will discuss later. 7 We emphasize that despite the repeated notation, the functions s00 and s0± controlling the Sommerfeld enhancement have numerically different values for the wino (Eq. (2.3)) and quintuplet (Eq. (2.21)).
Using these, we can evaluate, from which we conclude (2.24) With these modified forms for F 0 and F 1 , the LL result is then identical to that derived for the wino in Ref. [62]. For completeness, we restate it below.
where againα W = α W /4π and C A = 2. The first line in this result describes the two photon final state, arising from χχ → γγ, and is given in terms of a cross section that parameterizes the tree-level rate and a massive Sudakov logarithm, where s W = sin θ W , c W = cos θ W , and v = |v 1 −v 2 | is the relative velocity between the incident DM particles-notations we will use throughout. Substituting the Sommerfeld factors from Eq. (2.24) into Eq. (2.25), we see that the line cross-section is proportional to F 0 + F 1 = 4|4s 0±± + s 0± | 2 . This shows that the purely doubly-charged contribution to the line emission is a factor of sixteen larger than the purely singly charged, and also that the two contributions interfere. The equivalent result for the wino is 4|s wino 0± | 2 (again F 0 +F 1 using the wino equivalent values in Eq. (2.18)), which has the same form as the quintuplet when the doubly-charged contribution is turned off, although we caution that in the full result taking s 0±± → 0 does not reproduce the wino cross-section.
The second and third lines of Eq. (2.25) correspond to the endpoint, arising from χχ → γ +X, where the invariant mass of X is constrained to be near the lightcone. This contribution depends on an additional pair of logarithms and thresholds, associated with the jet (J) and soft (S) scales in the problem. In detail, (2.27) The extension to NLL proceeds identically. The expressions are more involved, but will be schematically identical to Eq. (2.25): all EFT functions will be identical between the wino and the quintuplet, with only the Sommerfeld contributions varying. To begin with, the differential NLL quintuplet cross-section can be written as 8 (2.28) The result is written in a similar form to the LL result of Eq. (2.25), with the exclusive twophoton final state on the first line, and the endpoint contribution on the last two. We note that Γ(x) is the Euler gamma function, not to be confused with the cusp anomalous dimensions introduced below. All functions in the result are identical to the NLL wino expression given in Ref. [64], except for F 0 and F 1 , which account for the various Sommerfeld channels. For completeness, let us first restate the elements common to the wino, again referring to Ref. [64] for additional details. Firstly, the evolution of the hard function is encapsulated in This result is written in terms of the first two perturbative orders of the β function and cusp 8 In this result we have set all EFT functions to their canonical scales. For instance, the weak coupling in the prefactor is evaluated at the hard matching scale µ 0 H = 2Mχ. A common technique for estimating the size of the theoretical uncertainty arising from neglecting higher order contributions is to vary these scales by a factor of ∼2. For this, the result with the scales unfixed is required, and can be obtained by extending the result in Eq. (2.28) to an arbitrary scale, exactly as done in Ref. [64]. anomalous dimension, as well as the ratio of the coupling between scales r H = α W (m W )/α W (2M χ ). The evolution of the jet and soft functions is contained in Here the ratio of scales are given by . Further, Θ J and Θ S are as defined in Eq. (2.27).
The last terms to be defined are those unique for the quintuplet. For those, we have and These expressions introduce r HS = r H /r S , as well as and further four functions Λ a−d . These functions are as follows, and are the same form as appears for the wino, with the functions ∆ written in terms of the polygamma function of order m, ψ (m) , as follows, (2.36) Finally, note that in the exclusive contribution to Eq. (2.28), we use the notation Λ → 1 to imply all of the functions in Eq. (2.35) are set to unity. Whilst the NLL expression is more involved than the LL result, the associated theoretical uncertainties are significantly reduced. This is demonstrated in the right of Fig. 1, where we show the cumulative, or integrated dσ/dz, taken from a given z cut to 1. The fact H.E.S.S. (or indeed any real imaging air Cherenkov telescope) does not have perfect energy resolution is represented by the fact the physically appropriate z cut is away from unity.
We will further explore these results in Sec. 5. Before doing so, however, we turn to the additional contribution to the spectrum that can result from bound states.

Bound State Formation
In this section we work out the rate of formation of relevant bound states, before considering the application of the SCET formalism to their annihilation in Sec. 4. The general formalism we employ is based on the methods of Ref. [70] for general non-Abelian gauge groups (see also Refs. [2,68]). Throughout this work, we consider only single-vector-boson emission in the dipole approximation. We first review the key equations and define our notation, then work out the form of the generators and potential with our basis conventions (Sec. 3.1), and use these results to evaluate the cross sections for bound-state formation via emission of photon (Sec. 3.2) and weak gauge bosons (Sec. 3.3). We note already that when SU(2) is broken, the velocity dependence of bound-state formation differs from that of Sommerfeld-enhanced direct annihilation; we will present a discussion of this issue and the uncertainties associated with the velocity distribution of the DM halo when we turn to our numerical results in Sec. 5.2.  If we label the different states in the multiplet containing DM as χ i where i runs from 1 to 5, then for capture of a χ i χ j initial state into a χ i ′ χ j ′ bound state with quantum numbers (nlm), where all particles have equal masses M χ and the emitted particle has color index a and can be approximated as massless, the amplitude for radiative capture into a bound state (stripping off the polarization vector for the outgoing gauge boson) is given by Ref. [70]: Let us define the various notations introduced in this amplitude. Firstly, α rad specifies the coupling associated with the radiated gauge boson: in our case, α rad = α for a radiated photon, 9 α rad = c 2 W α W for a radiated Z boson, and α rad = α W for a radiated W boson. T 1 and T 2 denote the generators of the representation associated with the χ i and χ j particles respectively, whilst f abc are the structure constants. (We emphasize that for the moment we are discussing the more general result where in principle χ i and χ j could be in different representations. Shortly, we will specialize to the case where T 1 = T 2 = T 5 appropriate for the quintuplet.) Finally, we have In the expression for J we have made the approximation of neglecting the momentum of the outgoing gauge boson. α NA is the coupling between the fermions and the t-channel gauge boson exchanged between them to support the potential, for the diagram where the bound state formation occurs through the emission of a gauge boson from the potential line. The two possible emission channels are depicted in Fig. 2. For example, when the bound-state formation occurs through emission of a photon or Z boson, the exchanged gauge boson will be a W boson, and so we will have α NA = α W . Theψ nlm andφ wavefunctions are the momentum-space wavefunctions of the final and initial states respectively (the corresponding real-space wavefunctions are labeled ψ nlm and ϕ).
The Y and J coefficients can be rewritten in position space as, 10 In subsequent equations we will suppress the r dependence of the position-space wavefunctions for notational convenience. We emphasize that these expressions tacitly assume the two particles are distinguishable; we will follow the conventions of Ref. [68] for the normalization of two-particle states, which can introduce a factor of √ 2 into terms involving transitions between states of identical and non-identical particles. 11 Breaking SU(2) can also introduce an additional multiplicative factor inside the integral for Y, arising from the propagators in the potential line from which the particle is emitted; for the case of photon or Z emission, these are W propagators, and the additional factor takes the form e −m W r [68]. We will work out the correct replacement in the case of W emission later in this section.
Throughout this section, when solving for the wavefunctions for the initial and final states given a specific potential, we will adopt the normalization conventions and numerical approach of Ref. [68]. (We note that there was a minus sign error in the equation for boundstate formation in the original published version of Ref. [68], which has since been corrected in an erratum.)

Generators and the potential
In general there is a degree of freedom to choose the basis for our generators, but since we are interested in transitions between two-body states whose constituents are mass eigenstates distinguished by their charges, it is convenient to use the basis discussed in App. A, where the χ 1 , χ 2 , χ 3 , χ 4 , χ 5 states correspond to states with electric charges +2, +1, 0, −1 and −2, as given in Eq. (A.9). It is important that the basis used to compute the bound-state formation rate and the basis used to compute the potential are identical; we will require the potential when solving for the initial-and final-state wavefunctions, which are relevant both for bound-state formation and for the Sommerfeld enhancement to direct annihilation. In this basis, we obtain for the generators: (3.4) 10 We follow the Fourier transformation conventions of Ref. [70]. 11 We emphasize that there is a subtlety here associated with the ordering of particles in the two-particle states, which can induce a sign flip that must be treated carefully. We discuss this in App. E.
The potential, up to terms corresponding to mass splittings between the two-body states (whose contribution is spelled out in App. A), can be written in the following form [4,72] Here N ij = 1 if i ̸ = j and 1/ √ 2 if i = j (this corresponds to the aforementioned change in normalization for two-body states composed of identical vs distinguishable particles), and the indices A, B run over {γ, Z, W + , W − }. The gauge couplings are included in the generators in this notation, following the conventions of Ref. [4]: explicitly, Note that the (−1) L+S factor arises from treating the ij and ji states as representatives of a single twobody state, using the conventions employed in method-2 of Ref. [72] and also discussed in App. E. (Here L denotes orbital angular momentum and S denotes spin; see Ref. [72] for a detailed discussion.) This sign was also discussed in the context of the potential for two-body states with net charge Q = ±1 in Ref. [2].
Thus, with the basis above, we obtain the following potential for the case with L + S even, where the 1st row/column corresponds to the χ ++ χ −− two-body state (χ 1 χ 5 ), the 2nd row/column corresponds to the χ + χ − state (χ 2 χ 4 ), and the 3rd row/column corresponds to the χ 0 χ 0 state (χ 3 χ 3 ): Note that the signs of the off-diagonal terms are opposite to the potential matrix given in Ref. [4] (and the analogue for the wino employed in Ref. [68]); this is a basis-dependent choice and either option is correct provided it is used self-consistently throughout the calculation. The effect of changing the basis in a way that modifies the signs in the off-diagonal terms of the potential is to flip the sign of one or more components of the resulting solution for the wavefunction; this compensates the changes in sign in generator elements in the new basis, when computing the bound-state wavefunctions. It is also possible to go beyond the tree-level potential of Eq. (3.6) and include NLO corrections. Especially in proximity to resonances, the resulting modifications to the Sommerfeld enhancement and bound-state formation rate can be substantial [5]. We employ the analytic fitting functions for the NLO potential calculated by Refs. [5,6]. Specifically, we make the following replacements in Eq. (3.6), where L ≡ ln(m W r): , m W r < 555/94 We will use the NLO potential for all calculations performed in this work. In particular, in addition to using it to compute the capture cross-sections in this section, we will also use it to compute the Sommerfeld enhancement appropriate for the direct annihilation discussed in Sec. 2. We note, however, that the relic abundance calculations in Refs. [2,3], which we rely on for our value of the thermal mass M χ = 13.6 ± 0.8 TeV, did not use the NLO potential. As we will see in Sec. 5, our findings as to the quintuplet will vary considerably across the thermal mass range, and so updating the calculation to NLO is an important improvement left to be done.

Bound state formation through emission of a photon
For the quintuplet it will often be possible to form a bound state through emission of W or Z gauge bosons, but let us first consider the case where the radiated particle is a photon, as this channel is available for all DM masses and provides the closest analogy to previous studies of the wino (e.g. Ref. [68]). In this case we have α rad = α, α NA = α W , and a = 3 (since the photon obtains its SU(2) couplings through the W 3 component). Let us also assume both incoming particles are in the same representation (as appropriate for Majorana fermions), and T denotes the generators of that representation, so we can write: From here, substituting in the generators above, we find the following non-zero matrix elements for bound-state formation which correspond to capture into the χ + χ − bound state component from the χ 0 χ 0 , χ + χ − , and χ ++ χ −− initial state components respectively. The equivalent matrix elements for capture into the χ ++ χ −− bound-state component are given by (3.10) Combining these matrix elements, and including a factor of √ 2 for the capture from the χ 0 χ 0 to χ + χ − state [68] to account for the differing normalization of states built from identical and distinguishable particles, we can write the cross section for bound-state formation as [70]: where k and ϵ(k) respectively denote the momentum and polarization of the outgoing photon; a CC subscript indicates the χ ++ χ −− component, C indicates χ + χ − , and N indicates χ 0 χ 0 . As previously, ψ and ϕ indicate the final bound-state and initial wavefunctions, respectively. Note that the potential of Eq. (3.6) is only accurate as written for two-body states with angular momentum quantum numbers L + S summing to an even value. If L + S is odd, the state is symmetric under particle exchange and cannot support a pair of identical fermions, and consequently the rows and columns corresponding to the χ 0 χ 0 state must be zeroed out. We are primarily interested in the behavior of quintuplet DM in the Milky Way halo, where on-shell charginos are not likely to be kinematically allowed (exciting the χ + χ − state requires 164 MeV of kinetic energy per particle, which for a Milky Way escape velocity of ∼500 km/s requires M χ ≳ 120 TeV). Consequently we will always assume the initial two-body state is χ 0 χ 0 (at large separation) and so has L + S even; this means the bound state formed by the leading-order vector-boson emission will have L + S odd (the dipole selection rule is ∆L = ±1, ∆S = 0). The appropriate potentials are used to compute the wavefunctions for the scattering state (Eq. (3.6) as written) and the bound state (Eq. (3.6) with the third row/column removed).

Bound state formation through emission of W and Z bosons
Since the SU(2) couplings of the Z boson are controlled by its W 3 component, we can re-use the expression for capture via photon emission in the case of the Z boson, with the replacement α → c 2 W α W in the prefactor, and with the momentum k now depending on the mass of the Z boson, Here the energy of the outgoing bound state is 2M χ + E n , where E n denotes the (negative) binding energy. 12 In particular, this process is forbidden when m Z exceeds the available energy (i.e. the kinetic energy of the incoming particles and the (absolute value of the) binding energy of the final state). The emission of W ± bosons is more complicated as it involves a different set of matrix elements and Feynman diagrams, corresponding to formation of bound states with unit charge. In particular, when a W boson is emitted from the potential, the t-channel propagator must be a mixed W − Z or W − γ propagator, which modifies the structure of the matrix element. Performing the Fourier transform of the mixed propagator, we find that the appropriate replacement (compared to the photon-emission case where the propagator involves only W bosons) is: where m 0 = m Z or 0 for the mixed W − Z and mixed W − γ propagator, respectively. Since the diagrams with these two propagator structures are identical except for the propagators and the coupling of the γ or Z to the fermion line, the sum of their contributions can be captured by inserting the propagator factor: (3.14) Now repeating the calculation from the photon case for the case where the emitted gauge boson is W 1 or W 2 instead of W 3 , and inserting the ζ(r) factors in the terms corresponding to W emission from the potential, we obtain the cross section for the Q = 1 case (the Q = −1 12 The convention here sets En = 0 as the rest mass energy of a pair of χ 0 s. For bound states that do not have a χ 0 χ 0 component, putting their constituents at infinite separation still leaves finite positive energy because of the charged/neutral mass-splitting, δ0. Thus, in a capture process, the bound state carries off energy Aδ0 − (|En| + Aδ0) = −|En|, where A is an integer that depends on the number of charged particles in the lightest component of the bound state. The second term, −(|En| + Aδ0) corresponds to the "ionization energy" necessary to separate the bound state into its constituents. We have neglected the subleading boundstate recoil kinetic energy. In this way we see that the boson emitted in the capture process has an energy independent of δ0. case is identical): Here the + + − subscript denotes the component in the χ ++ χ − state, and the +0 subscript denotes the component in the χ + χ 0 state, for the final Q = 1 bound state. Note that the phase-space factor k for the outgoing W boson must be modified to: As for Z emission, this process is forbidden when m W exceeds the kinetic energy of the incoming particles + the binding energy of the final state. The potential for the Q = ±1 sector, needed to derive the wavefunctions for the bound states, is similarly given by: where the first row/column corresponds to the χ ++ χ − state (Q = +1) or χ −− χ + state (Q = −1), and the second row/column corresponds to the χ + χ 0 state (Q = +1) or χ − χ 0 state (Q = −1). Again note that the off-diagonal terms disagree with Ref. [4] by a sign; this is due to our choice of basis.
In principle there may also be Q = ±2, 3, 4 bound states in the spectrum, which can be accessed by a series of transitions involving emission of W bosons. However, for the Q = 4 case, the only available state is χ ++ χ ++ (or χ −− χ −− in the Q = −4 case), and the potential is a repulsive Coulomb potential mediated by γ and Z exchange, which does not support bound states. For Q = ±3, the only available two-particle states are χ ++ χ + (χ −− χ − ), so again the potential is a scalar, and its value can be computed as We observe that for L + S even, this potential is always repulsive; for L + S odd, the potential vanishes in the unbroken limit and in the broken regime a residual repulsive potential remains.
In either case, we do not expect bound states (this analysis also accords with the discussion in Ref. [2]). The Q = ±2 case is more interesting. There are two relevant states: for Q = 2 they are χ + χ + and χ ++ χ 0 , with the former only being allowed for even L + S. The potential then reads as follows, for even L + S, where the 1st row/column corresponds to the χ + χ + state and the 2nd row/column corresponds to the χ ++ χ 0 state. This potential has an attractive eigenvalue that can support bound states, asymptoting to −3/r in the unbroken limit. For odd L + S only the χ ++ χ 0 state exists, which experiences no potential. Thus in addition to the Q = 0 ↔ Q = ±1 transitions through W emission already considered, the only bound-bound transitions we need to compute involving higher-charge states are Q = ±1 ↔ Q = ±2 (proceeding via W emission) and Q = ±1 ↔ Q = ±1, Q = ±2 ↔ Q = ±2 transitions via photon or Z emission.

Capture rate results
In Fig. 3, we show examples of the formation cross-section for bound states with different quantum numbers, corresponding to capture from various initial partial waves. Formation rates for Q = 1 bound states become non-zero at masses high enough that the binding energies (plus the kinetic energy of the collision) exceed the W boson mass. We observe that at most mass points, the dominant capture rate is to s-wave bound states, corresponding to the p-wave (L = 1) component of the initial state.
The overall size of the formation rate and its scaling with mass can be estimated analytically, as discussed in detail in App. F. To summarize, in the limit of high DM mass we expect  the leading rates for bound-state formation and direct annihilation to take the form: For M χ v ≲ m W , we expect the p → s capture cross section to experience a velocity suppression due to the p-wave initial state, which is parametrically of order (M χ v/m W ) 2 . In Fig. 4, we compare the dominant p → s bound-state formation rate for capture into the ground state with the inclusive direct annihilation rate; the latter is computed including the Sommerfeld enhancement but without any SCET corrections. We overplot the analytic estimates given in Eq. (3.20), with a p-wave correction factor of (M χ v/m W ) 2 for the estimate corresponding to the bound-state formation rate. We observe that at the thermal quintuplet mass (13.6 TeV), we expect the direct annihilation to dominate due to the p-wave suppression of the leading bound-state formation channel, but this suppression is lifted at high masses; furthermore, even at lower masses the bound-state capture rate may exceed the rate for direct annihilation at specific mass values (such as at the peak near 13.5 TeV). However, recall that these are inclusive rates; to understand the relative contributions to the line and endpoint spectrum, we must now understand how the bound states eventually annihilate to SM particles.

Bound State Annihilation
Having computed all the relevant bound-state formation rates, the second part of our calculation involves determining the differential branching ratios for bound states to decay, producing hard photons. To compute these we need the differential decay rate of the bound state to a final state including a photon, as well as its total decay rate to all SM particles. For the differential decay rate, we can recycle our EFT developed for direct annihilation as described in Sec. 2. The factorized form of the differential cross section remains identical to Eq. (2.14) which we reproduce here for convenience, To apply this expression to bound states, we will need to update the initial state wavefunctions encoded in F χ . For direct annihilation, F χ as given in the second line of Eq. (2.14), described an initial state of two free DM particles in the s-wave spin-singlet configuration. Here, however, our initial state is described by the two-body bound state wavefunctions computed in Sec. 3. The bound states can be classified according to their value of total orbital angular momentum L, total spin S, and charge Q, and we will need to track all bound states at a given mass. Beyond this, however, the differential cross section, dσ/dz, is again given by Eq. (2.5), and each of the associated objects such as the jet and soft functions are identical to those used in the direct annihilation computation of Sec. 2; this is the advantage of the EFT approach, the infrared (IR) physics is identical for direct and bound-state annihilation. To compute the decay rate, one needs to simply alter the details of the initial state, such as overall kinematic factors and a modified form of F χ . Nevertheless, as our interest is in the branching ratio of bound states to various decay rates, we will be computing ratios and will find the kinematic differences cancel (see Sec. 4.4), further increasing the similarity to the direct annihilation computation. We note, however, that to fully describe the possible end state of all bound states in the quintuplet spectrum, we would need to include additional operators in hard matching beyond the single operator we used for the direct annihilation given in Eq. (2.6). That operator described the annihilation of an L = S = Q = 0 initial state, so for the annihilation from states with L > 0, Q > 0, or S = 1, a new set of operators is required. Nevertheless, we will show that the contributions of the L > 0 and S = 1 states to the endpoint spectrum are suppressed, and the contributions from Q > 0 states can be captured within our existing framework, so that in fact the form of dσ/dz we have already computed is sufficient. Formalizing the logic above, the decay cross section into the hard photon can be written as where σ(χ 0 χ 0 → B + X us ) is the total production cross section for the L = 0, S = 0 state including any decays from shallower bound states and Γ B is the decay rate into all possible SM particles. We sketch the structure of the contributions to the endpoint spectrum from bound-state formation in Fig. 5, and work out the required ingredients in the rest of this section. We begin in Sec. 4.1 by understanding the general structure of the decay cascade that follows  . Schematic diagram that shows how two initial quintuplet particles evolve into various final states, as we work out in detail in this section. Endpoint photons are violet, whereas the soft photons that arise in the dipole transitions involving bound states are red. Capturing to quintuplet bound-states can give an additional source of endpoint photons. In this work, we include those that result from the "S = 0 Tower." In general, their contribution is suppressed compared to those from direct annihilation. As we discuss in Sec. 4.2.2, there is an enhanced capture rate to the "S = 1 Tower," including direct capture to the lowest-lying L = 0, S = 1 state, indicated by the thicker arrows. However, the decay to endpoint photons from this tower is power suppressed compared to S = 0 (hence the smaller endpoint photon on the right). This combination turns out to balance such that both towers give parametrically similar contributions. However, since this is a subleading overall component, only "S = 0 Tower" endpoint photons are included in our analysis.
capture into an excited state, arguing that excited states will typically decay (possibly via multiple steps) to an L = 0 state before annihilating to SM particles. In Sec. 4.2 we study the operators through which L = 0 bound states can annihilate to SM particles, and show that the contribution to the endpoint spectrum from S = 1 bound states is power-suppressed; thus our endpoint calculation focuses on annihilation from L = S = 0 states. In Sec. 4.3 we discuss how to compute the wavefunction factors needed to obtain the photon endpoint spectrum from decay (to SM particles) of a given L = S = 0 bound state. In Sec. 4.4 we compute the inclusive rate for decay via annihilation into SM particles for L = S = 0 states, while in Sec. 4.5 we describe how to calculate the rates for decay into lower-lying bound states, for bound states of arbitrary L, S. Key points of the calculation and several results are presented in Sec. 4.6. Ultimately, we employ these rates to compute the overall endpoint annihilation signal from decay of all L = S = 0 states to SM particles, taking into account the possibility to populate these states by decay from all shallower L = 0, 1, 2 states, as encapsulated in Eq. (4.2).

The decay cascade
Starting with an initial DM pair with a specific mass, we can expect capture into a number of metastable bound states characterized by their total orbital angular momentum L, spin S, and charge Q. These bound states have the option of either decaying into more tightly bound states with the same total spin but with |∆L| = 1 (in a single decay step), or annihilating directly into SM particles. 13 To compute the final annihilation spectrum into photons, we therefore need to know the branching ratios for various annihilation and decay channels. The decay of a shallow state into a low-lying "stable" state (by which we mean stable against decay to other bound states) may happen via several intermediate decay steps with their own branching ratios. We are therefore required to implement this cascade of decays to obtain the effective production cross section for a specific "stable" bound state, so we can then compute the signal from its subsequent annihilation into SM particles producing a hard photon. Determination of the full decay cascade requires three ingredients for the spectrum from bound states at a given DM mass: 1. the direct capture cross-section into all bound states; 2. the decay rate from each initial state to all more deeply-bound states; and 3. the rate for direct annihilation into SM particles. The first of these ingredients proceeds as discussed in Sec. 3. What remains to be computed is then the competition between the decay of one bound state to a deeper one, versus direct annihilation into the SM. 14 We will determine that in this section.
Before doing so, we can already provide an analytic estimate. The bare cross section for free electroweak DM particles to annihilate to the SM scales parametrically as σv ∝ v 2L α 2 W /M 2 χ , where L is the orbital angular momentum of the two-body initial state. The equivalent decay rate of a bound state to SM particles is related to this expression by replacing an incoming plane-wave wavefunction with the bound-state wavefunction. We can parametrically estimate this with two steps. Firstly, we replace v → α W as the characteristic momentum associated with the potential is p ∼ α W M χ and v = p/M χ in the nonrelativistic limit. Secondly, we must account for a multiplicative factor of (α W M χ ) 3 , which arises from the square of the bound-state wavefunction. (In more detail, as the wavefunctions are normalized by d 3 r |ψ| 2 = 1 and have support over the Bohr radius, a 0 , their characteristic value is |ψ| 2 ∼ 1/a 3 0 = (α W M χ ) 3 .) Thus we expect the decay rate of a bound state with orbital angular momentum L to SM particles to scale approximately as In contrast, a dipole-mediated decay to a lowerlying bound state scales as α α 4 W M χ independent of L; consequently, if such decays are allowed, they will generally dominate over annihilation to SM particles for L > 0. This argument suggests that unless dipole transitions to lower-lying states are forbidden, states with L > 0 will preferentially decay to L = 0 states before annihilating to the SM.
One might then ask whether the spectrum contains L > 0 states that have no allowed 13 Transitions where there is a change of spin or L by more than one unit are allowed, but suppressed, and we ignore them in this work; see e.g. Ref. [73].
14 Our discussion of this competition follows similar arguments in the literature, e.g. Refs. [68,74].
dipole transitions to more deeply bound states. For such states, the branching ratio for decay to SM particles via annihilation might indeed be important. However, we argue in App. F that this can only occur for very high-L states (beyond the range we consider in this work) for which the formation rate is likely to be negligible. We will thus assume that these states can be neglected, and restrict to L = 0 states when computing the endpoint photon spectrum from bound state decay. Similarly, in principle there can be stable Q = ±1, 2, L = S = 0 states in the spectrum. However, we find that the branching ratio to these states is very small (sub-percent) compared to the L = S = Q = 0 states, and so we neglect them in computing the endpoint photon spectrum. (In fact, at lower masses the charged bound state contribution will be exactly zero when either there is no charged bound state in the spectrum, or when those available cannot be accessed due to insufficient energy to produce an on-shell W .)

Leading power operators
In the direct-annihilation case, we expect s-wave annihilation (L = 0) to dominate. However, in order to support the annihilation of bound states with higher angular momentum we need operators that are suppressed by powers of the DM velocity. To see which will contribute, we consider the various structures that arise from a tree-level matching calculation. The tree-level amplitude for annihilation to a final state γ + X has contributions from s-, t-and u-channel diagrams, which give the following leading-order operator when expanded to O(v) As already mentioned, this operator is identical to that used for direct annihilation in Eq. (2.6) and supports a bound state with L = S = 0 (and both are given before a BPS field redefinition). As we will see, this is the primary operator required for computing the dominant contribution to the end-point spectrum from bound state annihilation. There is no operator at this order that supports S = 1, L = 0 bound state annihilation to gauge bosons at tree level. Such a state can annihilate to fermions and Higgs final states at tree level, but its contribution to end-point photons via bremsstrahlung is power suppressed in our EFT. These bound states however, can contribute substantially to the soft photon spectrum, as we will consider in Sec. 5.
For higher-L bound states with S = 0, there is a competition between decay to a state with lower L as compared with direct annihilation to SM particles. However, as discussed in Sec. 4.1, the decay to lower-L bound states always wins out for L > 0, so that only the decays of L = S = 0 states to SM particles remain relevant and the only operator we need is given in Eq. (4.3). Nevertheless, for completeness we provide the subleading operators in App. B. At the same time, there is no interference between the direct and bound state channels so that we may treat these cross sections separately. This is discussed in detail in App. C from an EFT perspective. If the widths of the bound states are parametrically much smaller than the separation in their energy, then we may also safely neglect any interference between the various bound state channels. This is essentially the narrow-width approximation. Accordingly, to determine the total spectrum it will suffice to sum over the cross sections for the direct channels and the allowed bound state channels individually.

Sub-leading power operators
As we saw in the previous section, there is no operator at leading power which supports an S = 1, L = 0 bound state annihilation into a hard photon. The only operators that support an S = 1 bound state are those which describe annihilation to fermions or scalars via an s-channel process. We can conceive of a hard photon emission from the final state SM Higgs or fermions; however, this is power suppressed by our SCET power counting parameter λ 1/2 (where λ = 1 − z) at the amplitude level. We can see this explicitly by looking at the emission amplitude of a hard (collinear) photon off a collinear fermion in the final state as in the diagram belowp p k The matrix element is then, where the first term is the tree level diagram and the second term comes from a single photon emission. By the power counting of SCET, we see that the collinear fields scale as A n,µ⊥ , f n ∼ λ, where λ is the expansion parameter of our EFT. The soft fermion field scales as f s ∼ λ 3/2 while the soft momentum p scales as n · p ∼ λ. This means that compared to the tree level, the hard photon emission is suppressed by a power λ 1/2 . If the fermion is ultra-soft, the suppression is enhanced to λ. Further discussion of this operator is provided in App. B. Including this operator in our analysis would be justified only if the production cross section for this channel compensates for the λ suppression (at the amplitude squared level), in order to be comparable to the S = 0 channel. Indeed this turns out to be true based on numerical calculations (see also App. F for an analytic estimate and discussion) and therefore, in principle we also need to include this sub-leading operator, in order to accurately compute the bound state contribution to the endpoint spectrum. However, a numerical analysis tells us that the leading S = 0 bound state channel, which will be the focus of the next subsection, is only a few percent of the direct annihilation cross section in terms of the contribution to the endpoint. The S = 1 bound-state contribution is power-suppressed relative to direct annihilation; it is only appreciable compared to the leading-power term from S = 0 bound states (whose formation rate is suppressed, by a different mechanism). So given the relative overall unimportance of the bound state contribution to the endpoint, we will not include the contribution from the S = 1 bound states here; a more precise calculation would need to account for this channel.

Wavefunction factors for bound state annihilation
In this section we compute the wavefunction factors F aba ′ b ′ relevant for the L = S = 0 bound states, which are needed to obtain the endpoint spectrum from the bound states' annihilation to SM particles. The definition of these factors remains identical to the case of direct annihilation given in Eq. (2.14), but now the operator is sandwiched between the L = S = 0 bound states.
As for the case of direct annihilation, the wavefunction factors will be evaluated at the IR scale of our EFT, which is the electroweak scale. At this scale, electroweak gauge symmetry is broken and the bound state wavefunctions are computed in terms of the broken eigenstates as in Eq. (2.21), but now for the bound state.
Once we have the bound state analogues of Eq. (2.21) in the broken basis, the next step is to relate these to the bound state wavefunction determined in Sec. 3. We start with the momentum-space representation of the bound state. In general for a two-particle bound state, we can express the bound state (which is an eigenfunction of the Hamiltonian) as where P is the momentum of the bound state, while 2k is the relative 3 momentum of the 2 particles making up the bound state. Further, ϕ P (k) is the momentum space bound state wavefunction whereas c † , d † are creation operators for the two constituents of the bound state.
The operators that we have in SCET are bilinear local operators with non-trivial Dirac structures. Based on the definition above we can evaluate the overlap of the bilinear operators with the bound state for the s-wave states. After SU(2) breaking, a general matrix element of a bound state will take the form (working in four-component notation for the moment and suppressing the color structure), Working in the bound state rest frame, we can expand this result to leading order in velocity, where ψ B (x) is the position space analogue of ϕ P (k), and we have introduced basis spinors η s and ξ r according to,v For a bound state in the spin singlet configuration, we can evaluate the basis spinors explicitly, and we are left with, For the computation at hand, we need to restore the color structure, which amounts to defining the bound state analogues of Eq. (2.21) which we used for direct annihilation. As our bound state is neutral, again there are only three objects to define where each of the bound state wavefunctions is the appropriate expression for that transition evaluated at the origin. We can then evaluate the wavefunction factors for the s-wave spin-0 bound state as follows,

Bound state decay rate into SM particles
In this subsection we compute rates for L = S = 0 bound states to decay into the SM. Because these states decay to gauge bosons, we can recast our previous results for photon emission through the L = S = 0 operator to obtain both the inclusive cross section and the differential branching ratio to photons. Taken together, these results will allow us to compute the hard photon spectrum from bound states, as dictated by Eq. (4.2). For the total (inclusive) decay rate for a given state, it suffices to look at the tree level cross section, since (as we explain next) there are no large logarithms induced due to loop corrections.
According to the KLN theorem [75,76], for non-abelian gauge theories, the cross section is IR finite only if we perform a sufficiently inclusive sum over both initial and final states. Since our initial states have a specific SU(2) color, one might expect to find IR divergences in the cross section when computing electroweak corrections. Since the electroweak symmetry is broken, these IR divergences should manifest themselves in the form of large logarithms ln(M χ /m W ). Most of the examples which demonstrate this violation of the KLN theorem in the literature consider light-like initial state particles as opposed to heavy time-like momenta considered in this paper. As we shall see, this is the key difference which removes the presence of large logarithms in inclusive cross sections for the case of heavy particle annihilation. Here, by heavy, we mean that the mass of the initial particle is of the same order as the hard scale in the EFT.
To demonstrate this explicitly, consider again the case of the wino but now imagine the wino to be a much lighter particle (mass ∼ electroweak scale) with a TeV scale (hard scale) energy. When we match the full theory onto an effective operator, the operator basis obtained after a tree level matching is fairly simple and reduces to (4.14) Note that in this case we have three soft functions in place of two since we can distinguish between the directions of the initial states. We are considering the inclusive case so that we sum over the colors of the final state gauge bosons. We can now look at the soft operators that we get at the amplitude squared level. As a single example, we can consider the interference term between Y 2 and Y 3 , which will contain Y ea n Y f b n Y f a n Y eb n . It is then clear that the soft Wilson lines do not cancel out due to the distinction between the n andn directions, if the initial state colors, a, b, are not summed over. On the other hand, if instead n =n = v, then the Wilson lines would cancel. This will render the soft function trivial and hence no Sudakov logs from KLN violation exist in this case. This implies that at NLL accuracy, for computing the inclusive cross section, we only need to consider the inclusive tree level cross section.
For the L = S = 0 operator, we need the differential branching ratio to the photon, in addition to the inclusive cross section. The differential decay rate is given by the following factorized formula which takes the same form as that of the differential cross section in Eq. (2.5) where A 0 is an overall kinematic factor. We do not need to know its explicit form since it will cancel out in the branching ratio.
To compute the inclusive cross section, we can make use of the stage 1 EFT given in Eq. (2.4) with all the functions evaluated to tree level where we have explicitly written out the convolution, and dΩ γ is an integral over the outgoing direction of the photon. This form is sufficient for computing the inclusive cross section, since we as explained earlier in this section, we do not have to resum any logs. An identical approach can be used to determine the inclusive decay rate of the L = S = 0 state, which at tree level decays purely to gauge bosons. Only several small alterations are required between the semi-inclusive (γ + X final state) versus inclusive cross section. We adjust the wavefunction factors contracted into the soft function, replace J γ by J n to allow for any final state gauge boson (we are no longer requiring a photon in the final state), and we remove the restriction on the phase space (as there is not observed endpoint photons) and therefore integrate over the full phase space. Further, we only need the various functions at tree level, and so we use For the inclusive case, we have a factor of 3 due to the color sum in the final state where we are now allowing for all gauge bosons instead of just a photon, although the function remains the same which is why we have retained the J γ notation. Then we are left with so now when we integrate over z and Ω γ we have The branching ratio, therefore can be obtained combining Eqs. (4.15) and (4.19) where the factor A 0 cancels out.

Bound state transitions
The previous subsections have established how to compute the photon spectrum from decay of L = S = 0 bound states via annihilation to SM particles. The remaining necessary ingredient is to determine how these states are populated through radiative capture and decays. The initial formation rate for bound states has already been discussed in Sec. 3; this subsection details the computation for shallowly-bound states to decay to lower-lying bound states.
Transitions between bound states, mediated by emission of a vector boson, can be computed using very similar expressions to those discussed in Sec. 3 for the initial bound-state formation. The are three salient differences: 1. the scattering-state wavefunction is now replaced with the bound-state wavefunction; 2. we must account for cases where the initial state has odd L+S and the final state has even L+S; and 3. we need address cases where the initial state has net total charge Q ̸ = 0. As discussed in Sec. 3, the second issue can be taken into account by modifying the potential used to compute the initial-and final-state wavefunctions.
The expression for the decay rate between states with total charge Q = 0 due to photon emission is given by a straightforward modification of Eq. (3.11), where as previously k is the momentum of the emitted photon and ϵ(k) is its polarization vector. Now ψ(r) is the wavefunction for whichever of the initial and final states has L + S odd, while ϕ(r) is the wavefunction for whichever state has L + S even (the dipole selection rule ensures that states connected by a dipole transition have opposite signs for L + S). We can make this simplification because the absolute value of the matrix element does not depend on the direction of the transition (although only one direction will have a positive k and hence non-zero available phase space). Consequently, once we have computed the matrix element for an L + S-even → L + S-odd transition (as in Eq. (3.11)), we can reuse the same matrix element for a transition in the reverse direction, only modifying the phase-space factors. The contribution from Z emission can be obtained by replacing α → α W c 2 W in the prefactor and replacing the momentum factor k as described in Eq. (3.12).
To see explicitly that the matrix element is invariant under time-reversal (up to conjugation), consider relabeling i ↔ i ′ , j ↔ j ′ in Eq. (3.1), and likewise swapping the notation for the initial-and final-state wavefunctions ψ nlm ↔ ϕ (that is, ϕ remains the wavefunction for the ij two-particle state, whether it is the initial or final state). The index relabeling has the effect of transposing the generator matrices; as the generators are Hermitian, this is equivalent to taking their complex conjugates. The relabeling of the wavefunctions applies complex conjugation to both J and Y, and additionally flips the sign of Y, as can be seen from Eq. (3.2). Consequently, we obtain: where Y and J are the functions computed from Eq. (3.2) for the ij → i ′ j ′ matrix element. Transitions between Q = 0 and Q = ±1 bound states are similarly related to the cross section for capture into Q = ±1 states, Eq. (3.15), with a rate given by: with ζ(r) as defined in Eq. (3.14), and phase space factor k as defined in Eq. (3.16). Here ϕ denotes the Q = 0 wavefunction and ψ denotes the Q = 1 wavefunction; in the case where the Q = 0 wavefunction has L + S odd, the terms containing ϕ N should be set to zero. For the Q = 1 sector there are no states containing two identical particles, so the set of allowed states is the same for L + S even or odd, although the wavefunctions for the two cases will differ due to the L + S-dependent potential.
Next, let us consider transitions mediated by photon or Z emission between two Q = ±1 states. In this case we obtain (suppressing the position dependence of the wavefunctions) for photon emission, and as above we obtain the Z-emission contribution by replacing α → α W c 2 W and modifying the phase-space factor k as defined in Eq. (3.12). In the expression above, initial and final states are distinguished by i and f subscripts, and this expression can be used for transitions where the initial state is either L + S-odd or L + S-even. Note that there is an explicit (L + S) i -dependent factor in the second line, which is necessary to ensure the consistency of the matrix element when the initial and final states are swapped. We discuss the details of its origin in Appendix E.
Finally, for transitions mediated by W emission between L + S-odd Q = ±1 states and L + S-even Q = ±2 states (as discussed above, we do not expect any L + S-odd Q = ±2 bound states), the decay rate is given by: where here ψ denotes the wavefunction for the Q = 1 state and ϕ denotes the wavefunction for the Q = 2 state, with subscripts labeling the components as above. The cross section is identical for transitions between Q = −1 and Q = −2 states, with the appropriate modifications to the wavefunction labels.

Key points from the bound state decay calculation
We now have all the ingredients for computing the differential decay rate of a bound state to a hard photon, so let us finally summarize our prescription (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 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.

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.

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σ/dz, both for the direct annihilation in Eq. (2.28), and also for the bound state contribution by combining the results of the previous sections with Eq. (1.1). For indirect detection, observables are sensitive to d⟨σv⟩/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, 15 This choice follows Ref. [62], and implies that our spectrum is normalized with respect to the line cross section, ⟨σv⟩ line ≡ ⟨σv⟩ γγ + 1 2 ⟨σv⟩ γZ , which is defined as the rate to produce two photons at exactly E = M χ . By construction, dN/dE will contain a contribution of exactly 2δ(E − M χ ) 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 ⟨σv⟩ i /⟨σv⟩ line , with ⟨σv⟩ i the cross section for that particular contribution. The rationale for anchoring our calculations to the line cross section is that χχ → γγ, which has a spectrum of exactly dN/dE = 2δ(E − M χ ), is a common experimental target, and therefore there are a wide number of existing constraints on ⟨σv⟩ line which we can then directly compare with. Further discussion of this point can be found in Ref. [62].
The spectra of line and endpoint photons produced by decay of the bound states is computed using the methods of Sec. 4. 16 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 15 Further discussion of the connection between spectra used in indirect detection and the corresponding field theoretic quantities can be found in, for instance, Ref. [77]. 16 In 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. [2,69].
(M χ v/m W ) 2L (although for sufficiently high masses, M χ ≳ 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 χ < 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 ⟨σv⟩ 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 χ . The figure makes clear a point already estimated earlier: direct annihilation generally dominates the production of line photons at E ≃ M χ 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 χ = 68.1 TeV, and therefore a reliable prediction at arbitrary masses must include this contribution.
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 χ = 13.6 TeV, and a mass where bound state contributions are significant, 68.1 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δ(E − M χ ) 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. [18]. In detail, we fix the relative width ∆E/E to 0.17 and 0.11 for E = 500 GeV and E = 10 TeV, respectively, and then vary logarithmically between these endpoints, freezing the ratio either side. From this we compute (dN/dE) 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. [62,64]: the endpoint contribution makes a considerable modification to the observed number of photons with E ∼ M χ , with the peak smeared spectra enhanced by 1.9 and 3.1 for M χ = 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 γZ and ZZ final states, and is a simple reweighting of the line cross section as For the W + W − final state, the tree level annihilation rate, with the Sommerfeld effect included can be computed as, As discussed for the case of the wino in Ref. [78], 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 ⟨σv⟩ W W /⟨σv⟩ line along with the Z contribution. In each case, to determine the spectrum of photons that result from these electroweak bosons we use PPPC4DMID [79] 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 γ ≪ M χ . 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 ⟨σv⟩ 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 → 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 → s capture rates. 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δ(E − M χ ) 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. [39], such rapid variations can lead to sharp features in the instrumental sensitivity to ⟨σv⟩ 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. [39]). The origin of this behavior seems to be interference between the different Sommerfeld factors, associated with the distinct mass eigenstates for the final annihilation: χ 0 χ 0 , χ + χ − , and χ ++ χ −− . 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 χ + χ − 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 χ 0 χ 0 and χ + χ − 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 dif- ference 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. [68]). We thus expect this behavior (sharp variations in the spectrum with mass) to be ubiquitous for larger SU(2) representations.

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 ≲ m W /M χ , which is true for halo velocities (v ∼ 10 −3 ) for M χ ≲ 80 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.
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 χ = 14 TeV, noticeable velocity dependence is present at v ≳ 4 × 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 ≲ 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 circ ≃ 240 km/s [80]. 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. [81]) but simulations including baryons have demonstrated the opposite behavior (e.g. Refs. [82][83][84]). 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. [85]). The escape velocity is determined to be ∼500 km/s at the location of the solar system [86][87][88]. 17 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. [89] 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 disp ∈ [130, 330] km/s. As a function of v disp , the magnitude of the relative WIMP velocity is drawn from the following 1D probability distribution, Here v disp is the RMS velocity for a single DM particle, which is equal to the three-dimensional velocity dispersion σ v,3d defined in Ref. [89] by Here f (r, v) is the speed distribution for a single DM particle at a distance r from the Galactic Center.
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.
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  Fig. 11. We see two different values of v disp plotted, plus a v fixed = 10 −3 whose capture rate we divide by.
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.

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. (5.1), the average DM-generated flux an instrument observes from a region of interest (ROI) of solid angle Ω ROI is, As defined, the flux has units of [counts/cm 2 /s/TeV/sr] (for a detailed discussion of the units, see Ref. [90]). 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 ρ χ 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. (5.6) 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. [62]. The approach in that work was to make use of the publicly available H.E.S.S. data in Ref. [18], where a search for χχ → γγ was center. As emphasized in Ref. [39], 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. [18] 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. [44]. Our rationale for adopting this mock analysis, however, is that Ref. [18] provided enough information that the full analysis they undertook can be performed reconstructed, as was shown in Ref. [62]. Later on we will provide a rough estimate of how sensitive more recent and upcoming analyses could be.
For the mock analysis, we fit the data provided in Ref. [18] to a combination of the flux in Eq. (5.6) and a parametric background model adopted by the experiment-full details are provided in Ref. [62]. If we have a prediction for ρ χ , this approach combined with our prediction for the spectrum can be used to obtain an estimated sensitivity for ⟨σv⟩ line . For this we adopt the Einasto profile [91] used by the H.E.S.S. analysis (based on Ref. [92]), with α = 0.17, r s = 20 kpc, and the normalization fixed to 0.39 GeV/cm 3 at the solar radius, r = 8.5 kpc. The resulting estimated constraint is shown in the left of Fig. 13. 18 We emphasize, that even though this is a constraint on ⟨σv⟩ 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. [18] 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. [62]). 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 ≃20%. We also note that there are several locations, such as just below M χ = 3 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. [39], 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 TeV mass range. Yet this statement is contingent on the form of ρ χ 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 ρ 2 χ flattening the profile has a marked impact on the flux, making this one of the dominant uncertainties in ρ χ 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": 19 We can then ask what choice of r c would raise the estimated constraint on ⟨σv⟩ 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 O(1 kpc) cores [94], however, measurements of stars in the Bulge seem to disfavor r c ≳ 2 kpc [95,96]. 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 [97]. 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 ⟨σv⟩ line ≃ 8.5×10 −27 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. [39] forecast a sensitivity at the thermal mass of ≃9.3 × 10 −28 cm 3 /s, almost a factor of ten better than used here. This is significantly more than the naive √ 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 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 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 [64] performed an identical mock analysis to ours for the NLL wino spectrum, estimating sensitivity at M χ = 13.6 TeV of ⟨σv⟩ line ≃ 8 × slightly stronger than the sensitivity to the quintuplet. Using the identical NLL spectrum, Ref. [44] then estimated that with 500 hours of data, CTA could reach ≃ 1 × 10 −28 cm 3 /s, a factor of eighty improvement. Assuming the same improvement for the quintuplet, CTA would be sensitive to ⟨σv⟩ line ≃ 1.1 × 10 −28 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 ρ χ out to almost 8.6 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 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.
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. [96]. 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 ± 0.8 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.

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 γ ∼ M χ 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 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 soonto-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 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.

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 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 χ = (χ 1 , . . . , χ 5 ) T , in terms of which the DM Lagrangian takes the following form 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 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 W , and is not a free parameter. Instead, M χ 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 χ = 13.6 ± 0.8 TeV [2,3]. 20 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 χ 0 , and two charged Dirac fermions χ + and χ ++ . At leading order, each of these states maintains a mass of M χ . However, radiative corrections to the charged states break this degeneracy, raising the masses of the charged fermions, singling out χ 0 as the lightest state and the DM candidate. These corrections have been computed, and in detail δ 0 = M χ + − M χ 0 ≃ 164 MeV, and δ + = M χ ++ − M χ + = 3δ 0 [1,98]. For most aspects of our calculations these mass splittings will be irrelevant and we will take δ 0 ≃ δ + ≃ 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δ 0 to the diagonal term in the potential matrix correspondingto the χ + χ − component of any state, 8δ 0 to the diagonal term corresponding to the χ ++ χ −− component, and similar, appropriate shifts for diagonal elements in the potential matrix corresponding to components of Q ̸ = 0 states (the shift is given by the difference between the rest mass of the state constituents and χ 0 χ 0 ).
We now turn to the interaction term in Eq. (A.1), 1 2 g Wχ T a 5 / W a χ. 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, The only part of Eq. (A.1) that remains undetermined is T a 5 , the three generators of SU(2) in the quintuplet representation. A convenient basis in which to specify T a 5 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. (A.2). In particular, as A µ 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 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 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 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 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 χ ij , and the quintuplet as χ ijkl , where the indices take values 1 and 2. Beginning with the wino, χ ij has three unique components, which we can embed into a vector as 21 We can use this representation to explicitly construct T a 3 as follows.
The key is that we know exactly how the generators act on χ ij as each index transforms in the fundamental. Accordingly, where T a F = σ a /2, with σ 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, κ) with κ ≪ 1, then the components of χ transform as From this, we can read off the action of an infinitesimal U on χ, and hence infer that we must have We can now identify χ 1 , χ 2 , and χ 3 as having charges +1, 0, and −1, yielding the expected spectrum in the broken phase. 22 The remaining components of T a 3 can be derived identically. This approach readily generalizes to the quintuplet. The five unique components of χ ijkl can be embedded into a vector as follows, To justify the charge assignments, we repeat the above analysis to find The √ 2 ensures that χ as it appears inχ / ∂χ/2 =χ ij / ∂χ ij /2 is canonically normalized. An identical argument explains the coefficients in Eq. (A.9). 22 The representation in the charge basis is not unique. For instance, the transformation χ 1,2 → e ±iϕ χ 1,2 leaves the charge assignments unchanged, but will introduce a phase into the off-diagonal W ± couplings. The same will be true for the off-diagonal quintuplet couplings.
Further, we can compute and T − 5 = (T + 5 ) T . In the main body we will exclusively work in the charged basis of Eq. (A.9), using the form of T a 5 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 χ ≃ 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 χ . For larger (smaller) values of M χ 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 χ ≤ 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 χ ≃ 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 σ ≃ (1.0 ± 0.3) × 10 −46 cm 2 [3] (see also Ref. [99,100]), beyond the reach of current searches at 13.6 TeV, and even next generation instruments such as LZ [101]. Nevertheless, the cross-section is a factor of ∼ 4 above the neutrino floor, and is in reach of generation-3 instruments such as DARWIN [102].
Detection at colliders is similarly challenging, but also potentially within reach of future instruments. Existing LHC searches reach masses of ∼ 270 GeV; even the high-luminosity dataset will only reach to ∼ 520 GeV [103]. 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 [104,105]. Given these two candidates have significantly lower thermal masses of 1 and 2.9 TeV respectively [4,[106][107][108], 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 √ s ≃ 35 TeV [3] (see also Refs. [109,110]). Taken together, in the short term indirect detection remains the most likely avenue for probing quintuplet DM.

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 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 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, 23 where it is understood that the subscript χ on v χ indicates that the velocity vector is the velocity of the state created/annihilated by the χ field (as opposed to theχ 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, 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. 24 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. (A.1) 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 γγ, γ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. 25 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 Yn Wilson line since we do not wish to consider sub-leading terms in the SCET power counting parameter.
We start again with our O(v 0 ) operator at the amplitude level before we dress it with soft Wilson linesχ where Γ is an arbitrary Lorentz structure. Consider the emission of an SCET ultra-soft gauge boson of momentum k off the initial χ orχ particle. First looking at the χ particle, p, p is the momentum of the χ particle. Let us now expand this result to O(v) (except inside the spinor). We see that apart from the usual Wilson line contribution, we also get two other relevant terms, If we sum the infinite series of gauge bosons with one of the propagators expanded out to O(v), we then have a structure where we have defined the following Hermitian operator, where P is the momentum label operator. If this term is included in a larger expression, the v · 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 25 One 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 [111], but this is not the case for the decay of L + S-odd states. In Ref. [68], one can see explicitly that for wino-onium, decays to W + W − are allowed for all combinations of initial-state L and S. emissions offχ to give us an effective operator, now dressed with soft Wilson lines We can once again use our Wilson line identity to write this in terms of our usual soft function We can then formally define a new object B s to explicitly separate out all the soft fields from the heavy DM fields where we have used Tr[T a T b ] = δ ab . Our operator now becomes In summary, we have an additional soft function which is explicitly suppressed in velocity. The Dirac structure for this operator is γ 0 γ 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, where 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 We can now move the |v| 2 factor to the DM wavefunction and instead redefine a soft operator The only term that contributes at one loop is the B since it is 0 at O(α 0 W ). So we can set all other Wilson lines to their tree level values. Explicitly, Thus, our operator at one loop becomes 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, where q + is the contribution to the final photon momentum from the soft function. The second term is proportional to M 2 χ and gives a power correction in M 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, Going to Laplace space and expanding out in the limit M χ → 0, we have 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 semiinclusive 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.

C Unstable Particle Effective Theory
In this section, we justify our use of Eq. (4.2) 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. [112]).
To begin with, if we have an intermediate resonance state we expect a propagator in our amplitude of the form where M * is a complex pole. If we write M * = M + iΓ/2, then we have the result We wish to work in a regime of narrow width, i.e. p 2 − M 2 ∼ ΓM ≪ M 2 , so that the propagator is well approximated by and we want to develop an effective theory with an expansion in the small parameter λ = Γ/M . For the case of inclusive decays of our DM bound state, ΓM χ ∼ α 5 W M 2 χ ; we are interested in the inclusive decay rate, so there are no large logarithms and the perturbative cross section begins at α 2 W , and due to the nontrivial wavefunction, we have an additional factor of at least α 3 W . The effective coupling thus scales as λ ∼ α 5 W ≪ 1. The hard scale in this process is just the resonance mass, which in our case is simply ∼M χ . We can now treat this as an HQET theory writing p = M χ 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 ∼ ΓM , we can immediately see that k ∼ Γ, so that we have a soft mode k µ ∼ (Γ, Γ, Γ) ≡ M χ (α 5 W , α 5 W , α 5 W ). Obviously, the question is how does this scale relate to the mass scale m W ∼ vM χ that we already have. Now the α W in the decay rate is evaluated at the scale M χ while the α W in the HQET scaling is at the scale m W ; however, we can see by the equations relating α W at the two scales that they are parametrically of the same order. If that is the case, then our mode has k µ ∼ M χ (α 5 W , α 5 W , α 5 W ) with k 2 ≪ m 2 W and hence can only be populated by a massless mode such as the photon.
Given the fact that in our case λ ∼ α 5 W , any corrections of the order λ 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 λ. Following [112], 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 µ only occurs at O(λ) and hence is severely suppressed. Additionally, since Γ ∼ α 5 W M χ , but the splitting between bound states, ∆E n ∼ α 2 W M χ , 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 ∼ Γ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 γ us indicating that a bound state was formed. In detail, the differential cross section is where M z is the measurement (in this case the photon energy) function on the final state particles, and 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, 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 ≫ ΓM ), then we can make the following substitution 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 Here B(p) represents a bound state with momentum p, and again this result is true at leading order in Γ/M , which forbids any communication between the production and decay states. If we then insert a factor of unity, 1 = d 4 p δ (4) (p − p f ), the result can be rearranged to yield, 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. (4.2).

D Proof of the Wilson Line Identity
In this section, we prove the identity involving soft Wilson lines used in Eq. (2.9) 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 where T is a generator in an arbitrary representation (we used 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 T a S v = S aa ′ v T a ′ , although this follows from the above by applying various inverses. In position space, the Wilson lines are defined as S v (x) = P e ig x −∞ ds v·As(vs) = P e ig v·x −∞ dv·y v·As(v·y) , where v · A s = v · A a s T a , with T a in the appropriate representation for the Wilson line, whilst P is path ordering andP indicates anti-path ordering. The variable s parametrizes the path along the light cone direction n from x to −∞. The statement in Eq. (D.1) is a generalization of the identity applied for QCD [113] for other group representations. The soft Wilson line obeys the equation Coming back to our question, let us define U a (x) = S v (x)T a S † v (x). We can then immediately see 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 ′ δ a ′ a . At the next order, (D.6) From here, we can see that the n th term will be Summing to all orders then proves our result.

E Subtle Signs in the Bound-state Formation and Decay Calculations
In the bulk of the paper, we have freely used results from Ref. [70] that are written in terms of two-body states of the form |ij⟩. However, in general our 2-body states for non-identical particles will in fact be combinations of the form 1 √ 2 (|ij⟩ + (−1) L+S |ji⟩). (This convention choice is also discussed in the context of Sommerfeld enhancement in Ref. [72], where these two approaches are labeled "method-1" and "method-2"; we largely adopt "method-2", where we treat |ij⟩ and |ji⟩ 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⟩ → |i ′ j ′ ⟩, we also need to include transitions between the |ji⟩ and |j ′ i ′ ⟩ states. In many cases this does not make a difference and it is adequate to represent states purely by one component |ij⟩. For example, if the |ij⟩ → |i ′ j ′ ⟩ and |ji⟩ → |j ′ i ′ ⟩ processes have equal rates, but |ij⟩ → |j ′ i ′ ⟩ and |ji⟩ → |i ′ j ′ ⟩ are forbidden (for example, this occurs if |ij⟩ = |0, ++⟩), then the combined rate is the same as what one would obtain from purely considering the |ij⟩ → |i ′ j ′ ⟩ 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⟩ → |+0⟩. Writing out the individual components of these states (labeled as 23 and 32, following the notation of Sec. 3), 26 the full matrix element should be: Now as calculated in Sec. 3, if ψ i and ψ f denote the initial-and final-state wavefunctions, we have: recalling that the "3" superscript is for γ, Z emission, depending on the value of α rad , the coupling of the emitted boson to the charged particle that radiated it. We also recall that the 26 One might be tempted to write the 23 state as |+ 0⟩ and 32 as |0+⟩. However, we reserve |+ 0⟩ for the full quantum state such that |+ 0⟩ = (|2 3⟩ + (−1) L+S |3 2⟩)/ √ 2.
terms with α NA correspond to emission off the virtual particles sourcing the potential. The α NA factor is thus the coupling of the virtual line to the WIMPs. For the case of capture by a radiated γ, α rad = α em and α NA = α W .
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.

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 boundstate 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. [68]), 27 The notation we employ follows Ref. [68]; η i and η f are potential eigenvectors which for the quintuplet are given by (in our basis) with corresponding attractive eigenvalues λ f = 5, λ 1 = 6, λ 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, λ 2 f α 2 W M χ /4. Lastly, theĈ 1 andĈ 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: In the wino and positronium cases studied in Ref. [68], there was only one eigenstate that experienced an attractive initial-state potential, and so the sum in Eq. (F.1) was trivial. There is a simple expression for |e πα W λ i /(2v) Γ(1 − iα W λ i /v)| 2 ≃ 2πα W λ i /v in the limit of small v (for positive λ i ), which scales purely as 1/v, and consequently in those cases (wino and positronium) σ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 scatteringstate wavefunction at the origin, and in the unbroken limit this wavefunction has the form i (I · η i )η i ϕ(λ i α W , r), where ϕ(α, r) is the solution to the scalar Schrodinger equation with an attractive Coulomb potential with coupling α. 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 χ 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. [68,114], 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 ≲ δ/m W , where δ is the mass splitting between the states; for δ = 164 MeV, we expect this behavior to hold roughly for v ≲ 2 × 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 W /M χ 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, λ i = 6 term contributing significantly, and with the coefficient I · η i replaced with δ 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 ≲ δ/m 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: Here we have employed the low-velocity approximation 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, In this unbroken limit, the effective branching ratio to the line (i.e. to γγ + half the branching ratio to γZ) should be given by (s 4 W + s 2 W c 2 W )/3 = s 2 W /3, so the line cross section should be: As for the bound-state formation, at higher velocities (v ≳ 2 × 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 approxi-mation holds (numerical calculations indicate the non-adiabatic cross section for bound-state capture, as estimated in Eq. F.1, can range between larger than the adiabatic result by a factor of ∼2 and smaller by a factor of ∼5, as v is varied). However, for v ≪ m W /M χ , we have the SU(2)-breaking effect that p-wave processes should be parametrically suppressed by a factor of order (vM χ /m W ) 2 , which at a 13.6 TeV mass can suppress the p-wave capture cross section by ∼2 orders of magnitude.
We can also study the cross-section for spin-singlet s → p capture to an n = 2, l = 1 state. In this case we have k = (25/16)α 2 W M χ , and σv = 2 13 παk 3 3 or in the adiabatic regime, 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 × 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 W /M χ , which renders it subdominant at the 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: Here the first row refers to the + + − state and the second to the + 0 state. The transition matrices are now,Ĉ In this case we must also replace the α prefactor in the cross section with α W . The attractive eigenvalue for the final state has λ f = 5, η f = {− √ 2, √ 3}/ √ 5. Again, the binding energy of the ground state is k = α 2 W λ 2 f M χ /4 = (25/4)α 2 W M χ (and since in the unbroken limit m W = 0, we do not need to include a kinematic suppression for the W mass). Then we obtain for the unbroken limit: This is for capture to the Q = +1 state; there is an equal rate for capture to the Q = −1 state. Note the α W prefactor (rather than α); 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 α W rather than α, for an overall capture rate (summing across all three channels) of ∼ 700π 2 α 3 W /(M 2 χ v), similar to the full s-wave direct annihilation rate. The primary difference between this p → s capture rate and the inclusive direct annihilation rate will arise from velocity suppression of the p → 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 χ , 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 /α 2 W M χ = −9/n 2 , −2.25/n 2 , and L-odd bound states with energies E n /α 2 W M χ = −6.25/n 2 . For spin-triplet states (S = 1) we expect L-odd states with energies E n /α 2 W M χ = −9/n 2 , −2.25/n 2 , and L-even states with energies E n /α 2 W M χ = −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 ′ = n, L ′ = 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 ′ = n, L ′ = 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 ′ = L − 1 will be available if 9/n 2 < 6.25/m 2 for some m with L ≤ m < n, i.e. if m < n 6.25/9 = n/1.2 is consistent with m ≥ L. This will be true for n > 1.2L, so the dangerous range is states with L < n ≤ 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 ≃ 0.14 in the Coulombic limit. The lowest-lying L = 6 spin-triplet state that is accessible via ∆L = 1, ∆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 spinsinglet 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 ′ = 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 ′ = L − 1 will be available except in the case where L < n ≤ 1.2L, which in the case of even L is potentially relevant for L ≥ 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 ≥ 5, L < n ≤ 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. [3]), 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 ≲ δ/m 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. [2] for general representations. They proceed by decomposing the two-particle state into eigenstates of isospin I (I = 1, 3, · · · , 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 λ = (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 λ = (N 2 − 1)/4. The isospin singlet corresponds to an L + S-even state with total charge Q = 0, and for the quintuplet λ = 6, as discussed above. In general, states with I < √ 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, ±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 λ i = (N 2 − 1)/4 and λ f = (N 2 − 5)/4. (Transitions amongst bound states can involve higher-I states; in particular I = 5 for the quintuplet, with λ = 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. [115]), 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. [2] in the simplified form: (F.12) We can now use tr(T a T b ) = T R δ ab , and also Thus, finally we obtain: where we will show how the C a1b J and C a1b τ 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 adj = 2, while d R = N (and in particular d adj = 3). Thus for SU(2) we obtain the coefficients: 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 λ. 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 ≤ N 2 −1 N 2 −5 L = 1 + 4 N 2 −5 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 isospinsinglet 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πα W λ i /v, in the low-velocity limit. This cross section is given for Majorana fermion DM by Ref. [2] as: Multiplying by the Sommerfeld factor gives: where in the final step we have assumed N ≫ 1.
Checking, for the wino and quintuplet this yields: For the quintuplet this agrees with our calculation above. This also agrees with the wino result from Ref. [68], accounting for our assumption that the adiabatic condition holds.

F.2.2 Capture to the ground state
At small velocities, from Ref. [2] we can read off the low-velocity bound state capture cross section to the n = 1 state as: This expression involves an average over initial states, with degrees of freedom denoted by g χ ; 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 χ = 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: where d adj = 3 counts the number of generators and arises from ab |δ ab | 2 = d adj .
In the limit of N ≫ 3, which is helpful for comparison against direct annihilation, we obtain the simplified result: (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. [2] in the low-velocity limit the s-wave (1st line) and d-wave (2nd line) contributions are: • s → p capture to the n = 2, L = 1 states collectively: direct annihilation rate ×0.06/N , in addition to any kinematic suppression of W/Z emission (up to a factor of 3α W /α), • d → p capture to the n = 2, L = 1 states collectively: direct annihilation rate ×1.8/N , in addition to any kinematic suppression of W/Z emission (up to a factor of 3α W /α) or velocity suppression due to the d-wave initial state (parametrically O(M 4 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 ≪ m W /M χ ) 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. [115]), 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 χ v/m 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 λ i /λ f > 1, leading to increasingly strong exponential suppression for largen states and avoiding the potentially unitarity-violating region of parameter space identified in Ref. [115].