Precise yield of high-energy photons from Higgsino dark matter annihilation

The impact of electroweak Sudakov logarithms on the endpoint of the photon spectrum for wino dark matter annihilation was studied intensively over the last several years. In this work, we extend these results to Higgsino dark matter χ10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\chi}_1^0 $$\end{document}. We achieve NLL′ resummation accuracy for narrow and intermediate spectral energy resolutions, of order mW2/mχ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {m}_W^2/{m}_{\chi } $$\end{document} and mW, respectively. This is the most accurate prediction to date for the yield of high-energy γ-rays from χ10χ10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\chi}_1^0{\chi}_1^0 $$\end{document}→ γ + X annihilation for the energy resolutions realized by current and next-generation telescopes. We also discuss for the first time the effect of power corrections in mW/mχ in this context and argue why they are not sizeable.


Introduction
In recent years much attention has been devoted to understanding the nature of the dark matter (DM) in the Universe. This type of matter is not accounted for in the Standard Model (SM) of particle physics and its existence is supported by compelling evidence from Galactic [1] to Cosmological scales [2]. Due to their accidental relationship with the electroweak scale, DM candidates known as weakly interacting massive particles (WIMPs) [3] have received most of the attention. Negative results from direct and indirect searches [4,5] of these WIMPs especially in the electroweak (EW) mass scale range suggest that more attention should be paid to less explored avenues such as the multi-TeV WIMP realizations.
A benchmark example for such heavy candidates is the Higgsino-like neutralino in supersymmetric extensions of the SM [6] and the limit of the pure Higgsino DM model, JHEP03(2020)030 which extends the SM by a single fermionic SU (2) doublet. Its mass has to be approximately 1 TeV if its cosmological abundance is to be explained by a freeze-out process. As a consequence of this large mass and the fact that in this scenario the Higgsino is then the lightest supersymmetric particle (LSP), the Higgsino DM model does not quite address the naturalness problem. The attractiveness of the model is its simplicity. However, direct searches for TeV-mass Higgsinos with the LHC are almost impossible and their scattering cross section with nucleons σ SI ∼ 10 −48 cm 2 [7,8] is below the so-called neutrino floor in direct DM detection experiments.
Nevertheless, DM models of the Higgsino type can be discovered indirectly if a line signature in the gamma-ray spectrum from, for example, the innermost region of our Galaxy, is observed. The quasi-monochromaticity of this part of the spectrum is a consequence of the kinematics of pair annihilation of non-relativistic WIMPs -a key prediction of many WIMP models including the Higgsino. This signal is particularly interesting because disentangling it from the uncertain astrophysical foregrounds is easier than for other type of spectra. Moreover, there is no known astrophysical mechanism that features all the aforementioned properties, rendering the observation of such spectral-line signals a smoking-gun discovery of annihilating DM.
Searches for this type of signature have already been performed by existing gamma-ray telescopes such as Fermi-LAT [9], H.E.S.S. [10] and MAGIC [11]. Particularly promising is the next-generation Cherenkov Telescope Array (CTA) [12], expected to improve the existing limits on spectral lines by one order of magnitude for TeV-scale gamma-ray energies. In these searches, though, it is assumed that the endpoint gamma-ray spectrum from DM annihilation is dominated by the fully-exclusive χ 0 1 χ 0 1 → γγ process. Under this assumption, the endpoint signal would be a perfectly monochromatic signal at the gamma-ray energy E γ = m χ as required by kinematics.
Even though the monochromatic approximation might be a reasonable assumption for some (light) WIMP models, it is not for generic models. The unavoidable finite energy resolution and the fact that only a single photon is observed at a time, imply that the proper observable is the flux that originates from the χ 0 1 χ 0 1 → γ + X process, where X is any type of (undetected) primary and secondary radiation that can be emitted in association with the observed gamma ray, and where the energy of the detected photon is close to the maximal value m χ within the energy resolution.
Theoretical predictions of spectral-line signals and the semi-inclusive photon spectrum near maximal photon energy for heavy WIMPs are not as straightforward as it might appear, since the perturbative expansion in the EW coupling constants breaks down. The effects giving rise to this are twofold. First, one has to take into account the so-called Sommerfeld effect, which is generated by the electroweak Yukawa force acting on the DM particles prior to their annihilation [13][14][15][16]. Secondly, for heavy DM annihilation into energetic particles, electroweak Sudakov (double) logarithms O((α 2 ln 2 (m χ /m W ))) are large and need to be resummed to all orders in the coupling constant [17][18][19][20][21][22][23]. The treatment of the Sommerfeld effect for Higgsino DM is well known. In this work, we focus on the resummation of large logarithmic corrections in the annihilation process χ 0 1 χ 0 1 → γ + X for Higgsino DM, using soft collinear effective theory (SCET).

JHEP03(2020)030
In recent papers we [22,23] and others [17][18][19][20][21] discussed these effects for the pure electroweak Majorana triplet model (wino-like neutralino) [24]. In those papers we established factorization formulas and obtained precise predictions for the photon yield for two different telescope energy-resolution regimes by using SCET methods [25][26][27][28]. We also provided the necessary functions for the evaluation of the factorization formulas at the next-to-leading logarithmic prime (NLL ) accuracy.
Here we extend these results to the case of the Higgsino DM, which is slightly more involved, due to the non-zero hypercharge of the Higgsino. As already mentioned, if the thermally produced Higgsino is to make up all the DM in the Universe, its mass should be approximately 1 TeV. This mass is indeed much larger than the masses of the EW gauge bosons and resummation of large Sudakov logs is necessary. Nevertheless, it is three times smaller than the mass of the thermally produced wino DM. As a consequence, m W /m χ mass corrections that are systematically neglected in the SCET leading-power resummations developed up to now might be as large as (or even larger than) the percent-level accuracy of NLL resummation. We thus include for the first time a quantitative discussion of the next-to-leading power mass correction for Higgsino DM (with direct applicability for wino DM as well).
The resummation of large Sudakov logarithms is achieved by breaking down the annihilation cross section into a few calculable factors representing the different scales and momentum configurations that build up the large logarithms. In order to prove such factorization formulas, some assumptions on the resolution E γ res of the photon energy measurement are necessary. In [23] we distinguished three parametrically different resolution regimes: 1 (1.1) Given the projected energy resolution of the upcoming CTA experiment [12] and the thermal Higgsino DM mass m χ ≈ 1 TeV, it is clear that the narrow and intermediate resolution regimes defined above are the most appropriate when computing the resummed cross section (see figure 1 in [23]). The article is organized as follows. In section 2 we briefly review the Higgsino model, summarize the SCET annihilation operator basis, and the factorization formula for the intermediate energy resolution. In section 3 we consider mass corrections. In section 4 we present the numerical calculation of the semi-inclusive annihilation cross section, after which we conclude in section 5. The main text is held short and focuses on non-technical results. In a series of appendices we give the potential used for the computation of the Sommerfeld effect, provide details of the factorization formula in the intermediate and narrow energy resolution regime and summarize the complete NLO expressions for all relevant functions that enter the factorization formula. 1 The theoretical framework for the "narrow" regime also applies to E γ res m 2 W /mχ including the classic γγ line with E γ res = 0.

JHEP03(2020)030
2 Endpoint spectrum of χ 0 1 χ 0 1 → γ + X The Higgsino model consists of the simple extension of the SM by an EW vector-like fermion SU(2) doublet with hypercharge Y = 1/2, such that one component is electromagnetically neutral after EW symmetry breaking (EWSB) [24,29]. The Lagrangian of the model is given by With our conventions (C.14) the choice of EW charges is such that the lower component of the multiplet is neutral, that is, in terms of components χ = (χ + , χ 0 D ), where the superscript denotes the electric charge. The mass eigenstates after EWSB are two selfconjugate (Majorana) particles (χ 0 1 and χ 0 2 ) defined in such way that χ 0 D = (χ 0 1 + iχ 0 2 )/ √ 2 and an electromagnetically charged Dirac (chargino) particle χ + .
A higher-dimension effective operator, for example L dim-5 = 1 Λ (χΦ)iγ 5 (Φ † χ) where Φ is the standard Higgs doublet, is necessary to provide the χ 0 2 particle with a slightly (at least O(100) keV) larger mass than the χ 0 1 particle. Otherwise, Z-boson mediated tree-level couplings of the Higgsinos to the light quarks would induce a large nucleon-DM cross section already ruled out by direct-detection experiments. On the other hand, the mass splitting δm = m χ − − m χ 0 1 between the charged and neutral component of the Higgsino doublet is induced radiatively after EWSB. At the one-loop order [30] Due to the Sommerfeld effect, the annihilation cross section is very sensitive to small variations of the mass splittings. However, the resummation of the Sudakov logarithms is insensitive to these variations. For instance, we checked that the effect of NLL resummation changes by less than 1% when the mass splitting of the neutral particles is varied by a factor of 10. Keeping this in mind, we will fix the chargino-to-Higgsino mass splitting to δm = 355 MeV, and further adopt δm N = m χ 0 2 − m χ 0 1 = 20 MeV for the mass splitting of the two neutral fermions. This value is small enough that dimension−5 operators responsible for it do not modify δm appreciably but, at the same time, is large enough that -for the range of DM masses considered here -χ 0 2 χ 0 2 cannot be produced for typical Galactic DM velocities v ∼ 10 −3 . We refer to [31] for a comprehensive discussion of the rich Higgsino mass splitting phenomenology.
The derivation of the factorization of the photon spectrum in the γ +X final state from DM annihilation near maximal photon energy 2 E γ = m χ has been thoroughly discussed in section 2.1 of [23], and holds for generic weakly interacting DM. For the sake of brevity, we will therefore provide here the model-specific expressions for the Higgsino and recommend the reader to consult [23] for theoretical background information.
Since the Higgsino multiplet has non-vanishing hypercharge, we must extend the basis (O 1-3 ) of short-distance annihilation operators with three more operators (O 4-6 ). For the annihilation process investigated here, only four out of the six operators are relevant (see JHEP03(2020)030 appendix C.1 for the complete set of operators) 3 where B µ is the SCET building block for the U(1) Y gauge field, and the spin-singlet matrix Γ µν is defined in [23]. In accordance with [23] the non-relativistic fields are represented in an unbroken-index notation. However, unlike in the wino case, in the unbroken Higgsino multiplet particles are not their own antiparticles. Thus, we adopt the nomenclature of [32] where the η v fields represent particles and ζ v the corresponding antiparticles. Additionally, for Higgsino DM. For j > 1/2 SU(2) multiplets the operators O 1 and O 2 are linearly independent. Hence, the factorization formula can be simplified by introducing the Wilson coefficientsC The photon energy spectrum in χ where S IJ captures the Sommerfeld effect, Γ IJ (E γ ) is the Sudakov-resummed annihilation rate, and the indices I, J run over Higgsino two-particle states χ i χ j . The Sommerfeld factor is computed with the leading-order potential in non-relativistic effective theory. The method of computation is discussed in [33] and the potential encapsulating the long-range force between the DM particles is given in appendix A. For gamma-ray energies E γ in a bin of (intermediate) energy resolution E γ res ∼ O(m W ) near the endpoint of the spectrum, 3 As in [32], the fields ζv and ηv are destruction operators for non-relativistic particles and antiparticles in the 2 and 2 * representations of SU(2), respectively (before EWSB). Note that after interchanging ζv and ηv, the operators remain the same (up to a sign in some cases). In contrast to [32] we do not include these redundant operators, and hence a factor of 2 appears relative to the operators defined in [22,23] for the wino case. With this convention the tree-level matching coefficient of the operator O2 is the same for Majorana and Dirac DM. 4 The overall factor of 2 is necessary in the method-2 computation of the Sommerfeld effect [33] for the annihilation of two identical particles to compensate for the method-2 factor 1/( √ 2) n id in (2.9). This factor of 2 has been missed in [22,23] and hence the absolute values of σv shown in figure 4 in [22] and figures 3-5 in [23] (published and arXiv version 1) must be multiplied by two.

JHEP03(2020)030
the Sudakov-resummed annihilation rate can be further factorized: This factorization formula resembles the one obtained in [23] and the meaning of the functions (and their arguments ω, µ and ν) is analogous:C are the short-distance coefficients of the hard annihilation processes, Z γ is the photon jet function describing the detected gamma-ray, J int is the jet function describing the unobserved hard-collinear radiation, and W describes the soft radiation. The indices W, Y = 3, 4 appearing in the soft functions and the photon-jet function are adjoint indices of the SU (2) . Specific expressions of these functions and their resummation are discussed in appendix C. Relative to the wino case, a new component due to the hypercharge of the Higgsino appears in (2.9).
In appendix B, we provide a more general version of the intermediate resolution factorization theorem, valid for general SU(2) multiplets, and more discussion of the functions appearing therein. We also provide the corresponding Higgsino-DM factorization formula for the narrow energy resolution regime E γ res ∼ m 2 W /m χ .

Mass corrections
The factorization formulas (2.8), (2.9) and (B.4), and the ones obtained in [22,23] are valid up to power corrections, i.e. terms proportional to v 2 , δm/m χ or λ ∼ m W /m χ and higher powers. The former two can safely be neglected as v 2 ∼ 10 −6 for DM annihilation in Milky Way-sized galaxies, and δm/m χ ∼ 10 −4 for the Higgsino and wino models. However, m W /m χ ∼ 0.1 × (1 TeV/m χ ) is larger than the percent-level accuracy of NLL resummed annihilation rates. In this section we shall investigate whether such linear O(λ), next-toleading-power mass corrections may be important. In order to assess the impact of power corrections on (2.8) we calculate the χ 0 1 χ 0 1 → γγ amplitude at the lowest, one-loop order in the coupling expansion. The computation of JHEP03(2020)030 χ 0 1 χ 0 1 → γZ would be similar. We find for the amplitude at vanishing relative velocity v = 0, expanded in λ = m W /m χ up to O(λ), the expression 5 In order to arrive at this result, we used feynrules [34] for the model implementation, FeynArts [35] for the Feynman rules and diagram generation, FormCalc [36] for amplitude processing, and PackageX [37] for the evaluation of the loop integrals. The finite and logarithmic pieces of the amplitude are all included in the Sudakov term Γ (11)(11) (Γ wino (00)(00) ). This can be verified by removing the m χ /m W and m W /m χ terms, squaring the resulting amplitude and comparing the result with eq. (313) in [23]. The term that is proportional to λ −1 = m χ /m W is the familiar one-loop Sommerfeld-enhancement factor.
Of particular interest for this discussion is the linear mass correction −2πm χ /m W × m 2 W /(24m 2 χ ). We find that this term is also associated with the non-relativistic dynamics of the problem. We verified this by showing that it arises exclusively from expanding the diagram displayed in figure 1 (and the one with the photon lines crossed) to subleading power in the potential loop momentum region. In the context of non-relativistic effective theory, the coefficient 1/24 originates from subleading-power potentials 6 (−1/8) and (at the squared-amplitude level) the matrix element of the dimension-8 derivative S-wave operator P( 1 S 0 ) introduced in [33,39] (+1/6). 7 Since in the non-relativistic theory leading-power contributions are O(λ −1 ), such O(λ) contributions should be counted as O(λ 2 ) corrections to the leading Sommerfeld enhancement. Together with the small coefficient 1/24, which results in a 3 × 10 −4 correction, we may conclude that we do not expect mass corrections to the leading-power factorization formula to degrade the percent level accuracy of the NLL resummation in the Higgsino mass range of interest.

Numerical results
As in [22,23] we turn our attention on the cumulative endpoint annihilation rate (4.1) 5 We also repeated this calculation for wino DM and obtained an almost identical result, iM χ 0 χ 0 →γγ For the numerical results given in this section we use the couplings at the scale m Z = 91.1876 GeV in the MS scheme as input: The upper panel in figure 2 shows the integrated spectrum σv (E γ res ) plotted as a function of the DM mass m χ , for the intermediate telescope resolution (2.8) set to E γ res = m W for definiteness. The displayed DM mass range includes the first Sommerfeld resonance. The different lines refer to the calculation of Γ IJ at tree level (black-dotted), the LL (magentadotted-dashed), the NLL (blue-dashed) and the NLL (red-solid) resummed expression for Γ IJ . The latter represents the result with the highest accuracy. For better visibility of the resummation effect the lower panel of the figure 2 shows the LL, NLL and NLL resummed annihilation rates normalized to the tree, i.e. Sommerfeld-only result. We see that resummation leads to a reduction of the cross section for large mass, as is generally expected for Sudakov resummation. In the low-mass regime around m χ = 1 TeV and below, however, resummation enhances the NLL and NLL annihilation rate.
This behaviour can be understood from the following observations. 1) At large masses the entries of the Sommerfeld matrix S IJ in (2.8) have similar magnitude and the sum over I, J is dominated by the (+−), (+−) term. The annihilation rate Γ (+−)(+−) in the diagonal χ + χ − → γγ channel starts at tree-level and has a standard series of exponentiated negative double-logarithmic corrections. On the other hand, for masses smaller than about 1 TeV the Sommerfeld effect is not very effective in mixing the various channels, and S (11)(11) ≈ 1, while the other elements are much smaller. Since, however, Γ (11)(11) is non-vanishing only from O(α 3 1,2 ), all terms in the sum over I, J can contribute equally to the annihilation rate and partial cancellations may occur. 2) Quite generally, at small masses it is not guaranteed that the leading logarithms dominate. Moreover, the leading logarithms in the neutral annihilation channels are positive, as we discuss below. This effect dominates over the negative interference term from Γ (11)(+−) and the Sudakov suppression of Γ (+−)(+−) , resulting in an enhanced annihilation rate at masses below 1 TeV.
The resummed predictions are shown with theoretical uncertainty bands computed from a parameter scan over the scales, with simultaneous variations of a factor of two of all scales (see also [23]). In the large-mass region the scale dependence decreases as higher orders are successively added from LL to NLL to NLL , but in the low-mass region the error band of the NLL prediction exceeds the error of the LL result and is very large in absolute terms. For small masses the inclusion of the non-logarithmic one-loop corrections to the hard, jet and soft functions is therefore necessary to gain control over the scale uncertainty, as is done in the NLL approximation. We then find that the resid-JHEP03(2020)030 ual theoretical uncertainty at the NLL order given by the width of the red-solid curve in figure 2 is small, and practically negligible for high masses. Numerically, for the two mass values m χ = 1 TeV (10 TeV The pattern of scale dependence at large masses is similar to the case of the wino [23], but the large NLL uncertainty at small masses was not seen there. An analysis of the analytic expressions for the resummed Higgsino annihilation matrix allows us to trace this behaviour to the following feature of the Higgsino model. The vanishing of the treelevel amplitude for the annihilation of a pair of neutral Higgsinos into γγ and γZ, which implies the vanishing of the tree-level annihilation matrix elements in the (00), (00) and off-diagonal (00), (+−) ((+−), (00)) components, 8 is caused by a cancellation between the short-distance coefficients of the three operators O 1,4,6 . This cancellation is not preserved by the electroweak renormalization group evolution, since the three operators have different anomalous dimensions and do not mix. 9 As a consequence there is a double logarithmic enhancement proportional to L 2 = ln 2 (4m 2 χ /m 2 W ) in the one-loop amplitudes despite the absence of a tree amplitude. Then the lowest non-vanishing order in Γ NLL (00),(00) , which is O(α 4 1,2 ), 10 already carries a L 4 enhancement. This is in stark contrast to the wino model, where Γ (00),(00) contains at most L 2 at O(α 4 2 ) (see appendix E of [23]). The different logarithmic structure is related to the hypercharge of the Higgsino, which allows the decay The higher powers of logarithms in the neutral annihilation matrix channel, which is relevant at low masses, is also the origin of the large scale dependence of the NLL approximation. Defining l µs ≡ ln(µ 2 s /m 2 W ) and l µ h ≡ ln(µ 2 h /4m 2 χ ), we have for the first non-vanishing, O(α 4 2 ) contribution Γ NLL (00),(00) =α where we write explicitly only the terms relevant to the discussion. The existence of a L 4 term implies that the coefficients of L 2 and L are dependent on the matching scales, consistent with the fact that these coefficients are not yet properly summed in the NLL approximation. We then find that the large scale dependence is caused by the π 2 -enhanced terms shown in (4.2), which in turn stem from the imaginary parts of the one-loop anomalous dimensions. The variation of these single-logarithmic terms under scale variation is larger than the scale-independent part of Γ NLL (00),(00) , which causes an O(1) scale dependence at small masses. 11 Adding the one-loop non-logarithmic terms to the hard, jet and soft JHEP03(2020)030 functions at NLL removes the scale-dependent terms shown in (4.2), which explains the dramatic reduction of the theoretical uncertainty when going from NLL to NLL . The same large scale-dependent terms are also present in the LL approximation. However, in this case there is an accidental cancellation of large scale dependence between the (00), (00) and (00), (+−) contributions to the sum in (2.8), resulting in an accidentally small, and highly asymmetric scale uncertainty, as seen in figure 2. None of these observations is relevant to the high-mass regime, where the sum is by far dominated by the (+−), (+−) component, which has a very small scale dependence already at NLL. Neither are they for the wino model, where the large scale dependent terms do not appear. Furthermore, the high-mass regime sets in earlier for the wino, because the Sommerfeld potential is stronger than for the Higgsino.
As for the case of wino DM, it is instructive to investigate whether the narrow (E γ res ∼ m 2 W /m χ ) and the intermediate (E γ res ∼ m W ) resolution results can be matched to provide an accurate prediction for the entire range from E γ res ∼ 0 to E γ res ≈ 4m W . In figure 3, we show the annihilation cross sections for the narrow (blue-dotted) and the intermediate We observe that there is a wide range of E γ res , where the two resolutions match very closely. This implies an accurate theoretical prediction for the photon energy spectrum in JHEP03(2020)030 DM annihilation in the entire range from E γ res ∼ 0 to E γ res ≈ 4m W . A similar matching was performed in the wino DM case and we found the same degree of matching for the two resolution regimes. An in depth investigation of the structure of the large logarithms was performed for the case of wino DM in [23], which explained the high degree of matching of the two resolution cases. A similar structure holds for the Higgsino.
There is a steep rise in the narrow resolution cross section, which occurs at E γ res m 2 Z /(4m χ ). Above this value, the γZ contribution cannot be resolved, resulting in the sharp increase of the semi-inclusive rate. Since the unobserved jet function for the intermediate resolution regime is computed assuming massless particles, it misses this effect. The invariant mass of the narrow resolution unobserved jet function also passes through the W + W − , ZH and tt thresholds, which are however invisible on the scale of the plot.
Due to the high-precision agreement between the two resolution regimes over a wide range of E γ res , we can conclude that in figure 2 the integrated photon energy spectrum of the narrow resolution case would look indistinguishable from the intermediate resolution results, provided the chosen value of E γ res lies somewhat above m 2 Z /(4m χ ) and below 4m 2 W /m χ . After performing a similar parameter scan as in the intermediate resolution case to determine the theoretical uncertainty, we find that the reduction of the scale uncertainty with increasing accuracy is comparable for the narrow resolution case.

Conclusions
The differential annihilation cross section for Higgsino DM annihilation into γ + X near the endpoint was calculated in this paper with the greatest accuracy to date (O(2%) at m χ = 1 TeV and less at higher masses). This computation systematically includes the resummation of the Sommerfeld effect and electroweak Sudakov logarithms at the NLL order and is valid for intermediate and narrow energy resolutions. To achieve the above-mentioned accuracy, it is necessary to perform resummation at NLL accuracy. Our predictions are appropriate for dedicated searches for Higgsino DM (within the thermal production hypothesis) using next-generation gamma-ray telescopes. Electroweak Sudakov resummation decreases the rate of photons by approximately 3% for a Higgsino mass of 1 TeV. For larger and smaller masses resummation leads to suppression (45% at 10 TeV) or enhancement, respectively (numbers refer to E γ res = m W ). We also discussed for the first time in the context of indirect DM detection the impact of power-suppressed mass corrections on the endpoint spectrum. In the case of wino and Higgsino DM we find these corrections to the one-loop amplitude to be of O(m 2 W /m 2 χ ) relative to the leading Sommerfeld effect and thus compatible with the accuracy achieved by NLL computations. The largest theoretical uncertainty for the spectrum is expected to be associated with NLO corrections to the Sommerfeld potential, which has recently been computed for the wino [38] but is not yet known for Higgsino DM.
The Higgsino annihilation cross sections σv (E γ res ) ∼ 2×10 −28 cm 3 /s (for m χ = 1 TeV) predicted here are below the limits obtained by existing gamma-ray observatories but will be probed by the CTA experiment. Concretely, the H.E.S.S. experiment quotes an upper limit (2σ) on σv γγ of 4 × 10 −28 cm 3 /s [10] for their most aggressive assumptions on JHEP03(2020)030 the DM mass distribution of the Milky Way, which translates to about σv H.E.S.S. limit γ+X ≈ 8×10 −28 cm 3 /s for the semi-inclusive rate with energy resolution E γ res ≈ 100 GeV. With the help of our results, future measurements of CTA can be converted into precise constraints on the Higgsino mass given the dark matter profile of the galaxy.

Acknowledgments
We thank Clara Peset for very useful discussions. This work was supported in part by the DFG Collaborative Research Centre "Neutrinos and Dark Matter in Astro-and Particle Physics" (SFB 1258).

A Sommerfeld potential
In our numerical evaluations, the method-2 Sommerfeld matrix S IJ (see [33]) is obtained by solving the Schrödinger equation with the spin-singlet potential The indices are ordered in the following way: (11), (22) and (+−). We added the contribution of the mass-splitting matrix δm IJ . As mentioned in the main text, there is no interaction between the above three two-particle states and the mixed (12) state.

B.1 Intermediate resolution regime
The derivation of the factorization theorem is independent of whether the DM particle is charged under hypercharge or not. Using the result from [23], we can thus write the general form of the Sudakov-resummed annihilation rate for the intermediate resolution case as As already discussed in section 2, we now augmented the operator basis to include hard annihilation into the hypercharge gauge boson. The anti-collinear function Z Y W γ (photon jet function), the hard-collinear function J XV int (unobserved-jet collinear function) and the soft function W ij IJ,V W XY are generalizations of the corresponding jet and soft functions defined in [23]. In particular, the SU(2) ⊗ U(1) Y indices W and Y can now adopt the values 3 and 4, where 3 refers to the SU(2) gauge boson W 3 and 4 to the U(1) Y gauge JHEP03(2020)030 boson B. It is convenient to split the unobserved jet function into an SU(2) and a U(1) Y component, as follows: was already obtained in [23], while J

U(1) int
is new to the case of Higgsino DM. The results are presented in appendix C.3. With (B.2) we can reduce the number of relevant indices of the soft function by introducing Here X, V are summed from 1 to 4. For the pure Higgsino with isospin j = 1/2, we can also make use of the degeneracy of operators O 1 and O 2 in (2.6), from which the factorized Sudakov annihilation rate (2.9) immediately follows.

B.2 Narrow resolution regime
Assuming the energy resolution E γ res ∼ m 2 W /m χ puts us into the narrow resolution regime. The hierarchy of scales then changes to E γ res m W , m X ∼ m W , which means that the unobserved jet now has collinear rather than hard-collinear virtuality. Furthermore, real soft gauge boson radiation is power suppressed making it convenient to write the soft function as soft Wilson coefficients D. A detailed discussion of the factorization theorem for wino DM in the narrow resolution case is provided in [22] and [23]. Since there are no conceptual differences for Higgsino DM, we will not repeat it here. The Sudakov-resummed annihilation rate Γ IJ for the narrow resolution case is given by (the Sommerfeld factor S IJ remains the same) where the SU(2) ⊗ U(1) Y indices summed over are V, W, X, Y = 3, 4.

C NLO expressions
In this appendix we provide all expressions that are required for the evaluation of the factorization formulas (2.9) and (B.4) at the one-loop order as is needed for the NLL resummation. We also provide the necessary results for resumming the various functions.
Since there are no conceptual differences with respect to the resummation in the wino DM case, we refer to [23] for an in depth discussion of the systematics of resummation and will keep the present discussion rather brief. Also, all function definitions have already been presented in [22] and [23] and will not be repeated here.

C.1 Operator basis and Wilson coefficients
In the general case where the field χ in (2.1) is a (2j + 1)-plet of SU(2) with non-vanishing hypercharge, the operators allowed by the symmetries of the SM are The derivation of this operator basis is as in [23]. As argued there, operator O 3 is irrelevant for the χ 0 1 χ 0 1 → γ + X process. We also find that O 5 is irrelevant as it is a spin-triplet operator similar to O 3 . The spin-triplet operators do not contribute, since there is no spintriplet χ 0 1 χ 0 1 initial state, and the Sommerfeld-enhanced scattering prior to annihilation does not change the spin.
The one-loop short-distance coefficients of operators O 1 and O 2 have already been obtained in [23]. For completeness we provide the Wilson coefficients for the entire operator basis (C.1)-(C.6). The Wilson coefficients are computed from matching the full-theory amplitude to the effective theory amplitude. The Feynman diagrams are calculated in the unbroken SU(2) ⊗ U(1) Y gauge theory. For general j and Y , the coefficients are given by

JHEP03(2020)030
Here c 2 (j) = j(j + 1) is the SU(2) Casimir of the isospin-j representation and n G = 3 is the number of fermion generations. 12 Note that C 3 (µ)| Dirac in (C.7) is specific to a Dirac SU(2) ⊗ U(1) Y multiplet. If, on the other hand, we assume a Majorana multiplet, we find Let us comment on a subtlety in the computation of C 3 . Naively, one would expect that there should be no counterterm contributions as the tree-level contribution cancels. However, as also observed for the corresponding quarkonium calculation in QCD [40], the vanishing tree-level result stems from a cancellation between an s-channel diagram and the t/u-channel diagrams. Since DM mass counterterm insertions exist only for the tand u-channel diagram, a mass renormalization contribution survives, which is required to obtain a finite result for C 3 . One might wonder how this result can be obtained in bare perturbation theory, as the tree-level result vanishes, and hence there is apparently no bare mass to substitute by the renormalized one. In bare perturbation theory, the appearance of the counterterm is related to a subtlety in the matching procedure. When performing one-loop on-shell matching we must set the mass of the external particles to their renormalized on-shell values at the corresponding loop order. The tree-level diagrams then contain the bare mass from the explicit mass in the propagator and the renormalized mass from momenta after applying on-shell kinematics. The above mentioned cancellation then leaves over the one-loop difference between bare and renormalized mass in the t-and u-channel tree diagrams and the same result as before is recovered in bare perturbation theory.
The operator O 5 is irrelevant for χ 0 1 χ 0 1 annihilation, but we also find its coefficient to be zero at the one-loop order. This is a consequence of the Landau-Yang theorem. Although it does not hold in non-abelian gauge theories (see, for instance, [40,41]), the violation arises due to the fact that the final state bosons carry an internal quantum number and structures can be constructed involving the group structure constant. At least to the oneloop order, there are no such Feynman diagram structures for O 5 as one of the final-state gauge bosons is abelian.
Next we discuss the renormalization-group (RG) evolution of the short-distance coefficients. Recall that we are considering the annihilation process χ 0 1 χ 0 1 → γ + X for Higgsino DM, for which the operators O 3 and O 5 are irrelevant. Hence, we will disregard C 3 (µ) and C 5 (µ) from now on. The evolution of the vectorC = (C 1 ,C 4 ,C 6 ) defined in (2.6) and (2.7) is given byC The matching coefficient for O3 was previously given for j = 1 and Y = 0, i.e. the pure wino, in [20].
The result given there differs from ours. We attribute this difference to an opposite sign in one of the diagrams in [20] -namely T 5b , and missing mass renormalization (counterterm) diagrams. However, as mentioned before, O3 does not contribute to the annihilation rate into photons.

JHEP03(2020)030
where n i,SU (2) and n i,U (1) give the number of SU(2) and U(1) gauge fields in operator O i . Eq. (C.10) is solved numerically to obtain the resummed short-distance coefficients (C.9). Most anomalous dimensions in (C.11) have already been given in [22] and will not be repeated here. The new ones relevant for the evolution ofC 4 (µ) andC 6 (µ) are γ U(1) (α 1 ) = γ The result for the photon jet function Z 33 γ was already given in [22,23] and will not be repeated here. For the case of Higgsino DM, which has non-vanishing hypercharge, we also need to take into account the index combinations Z 34 γ , Z 43 γ and Z 44 γ , which can be obtained by adapting Z 33 γ accordingly. To do so, it is helpful to remark that only diagram (a) ((b)) contributes to Z 34 γ (Z 43 γ ), while Z 44 γ does not receive contributions from (a) and (b). Hence, the Wilson line part of Z 34 γ and Z 43 γ is multiplied with a factor of 1/2 with respect to Z 33 γ while Z 44 γ does not have a Wilson line contribution at all. To make the computation more transparent, it is convenient to write down the rotation from the weak basis to the mass JHEP03(2020)030 of the gauge fields. One can now use (C.14) to compute the remaining photon jet functions Z 34 γ , Z 43 γ and Z 44 γ by carefully analyzing the position of the cut in diagrams (a), (b) and (c) and by keeping track of whether the γ or the Z boson originated from a W 3 or a B boson in the unbroken theory. We find The discussion of the RG and rapidity RG evolution will be limited to Z 34 γ = Z 43 γ and Z 44 γ . The resummation of Z 33 γ was discussed in [23] in great detail, which allows us to keep the following analysis rather brief. The RG equations are given by with the anomalous dimensions The anomalous dimensions γ

JHEP03(2020)030
Since Z 44 γ is independent of the rapidity scale ν, we only need to consider the rapidity RG equation for Z 34 γ . It reads where the rapidity anomalous dimension γ ν Z 34 γ is given by Solving (C. 18) and (C.22), we can compute the resummed photon jet functions where the integrals in the exponents are computed numerically.

C.3 Unobserved jet function C.3.1 Intermediate resolution
In addition to the SU(2) gauge-boson jet function J SU(2) int (p 2 , µ) obtained at NLO in [23], the factorization formula (2.9) also requires its U(1) Y counterpart. It is given by where l = 1/(e γ E τ 2 ) and the explicit result reads The corresponding RG equation is the ordinary differential equation with the Laplace-space anomalous dimension is needed at the one-loop order for NLL resummation: The RG equation (C.27) is solved by where µ j ∼ 2m χ m W is the natural scale of the hard-collinear jet function. To return to momentum space, we perform a standard inverse Laplace transform (remembering that τ 2 = 1/(e γ E l)). We find for the resummed unobserved jet function J where U (µ j , µ) is taken from (C.30) and J U (1) int (p 2 , µ j ) on the right-hand side of (C.31) is given by (C.24). The integral in the evolution factor U (µ j , µ) is evaluated numerically.

C.3.2 Narrow resolution
As for the photon jet function discussed in section C.2, the unobserved jet function in the narrow resolution case receives contributions from the index combinations J 33 (p 2 ), J 34 (p 2 ), J 43 (p 2 ) and J 44 (p 2 ). The function J 33 (p 2 ) was already computed in [22] and a detailed discussion of its derivation is presented in [23]. The result is a complicated function of the masses of the SM particles, the invariant mass of the jet and of the virtuality and rapidity scales µ and ν, respectively. The computation of J 34 (p 2 ), J 43 (p 2 ) and J 44 (p 2 ) is similar to the case of the photon jet function. One has to be careful which functions receive contributions from Wilson line diagrams and whether the gauge boson crossing the cut originates from a W 3 or a B boson, in order to determine the correct prefactor using (C.14). Since J 34 (p 2 ) = J 43 (p 2 ) we will content ourselves with presenting the results for J 34 (p 2 ) and J 44 (p 2 ). They read for J 34 (p 2 ) and (ω) = δ (ω) , (C.46) where we introduced n ij IJ = (−1) δ I(00) δ i4 (−1) δ J(00) δ j4 with (00) = (11) or (22), (C.47) to allow for a more compact notation. The resummation of the soft function is conceptually analogous to the resummation of the soft function in the case of wino DM [23], to which we refer to for a detailed discussion. The Laplace transform of the soft function and its inverse are defined as (ω, µ s , ν) . (C.64) As in [23], we did not include the rapidity evolution factor in (C.64), since we evolve the photon jet function in ν from ν h to ν s , which makes the soft function rapidity evolution factors unity (see (C.58)).

C.4.2 Narrow resolution
In this appendix, we collect the soft functions appearing in the factorization theorem for the narrow resolution case (B.4). The definition and computation of these terms is analogous to the wino DM case and we refer the reader to [22] for more details. The expressions take the following form: