Resummed photon spectrum from dark matter annihilation for intermediate and narrow energy resolution

The annihilation cross section of weakly interacting TeV scale dark matter particles $\chi^0$ into photons is affected by large quantum corrections due to electroweak Sudakov logarithms and the Sommerfeld effect. We extend our previous work on the resummation of the semi-inclusive photon energy spectrum in $\chi^0\chi^0\to \gamma+X$ in the vicinity of the maximal photon energy $E_\gamma = m_\chi$ with NLL' accuracy from the case of narrow photon energy resolution $E^\gamma_{\rm res}$ of order $m_W^2/m_\chi$ to intermediate resolution of order $E^\gamma_{\rm res} \sim m_W$. We also provide details on the previous narrow resolution calculation. The two calculations, performed in different effective field theory set-ups for the wino dark matter model, are then shown to match well, providing an accurate representation up to energy resolutions of about 300 GeV.


Introduction
High-energy photons may constitute an important signal for the particle nature of dark matter (DM) through the pair annihilation of DM particles. In order to distinguish the DM component from the astrophysical γ-ray background, one searches for the line signal of the two-body annihilation χ 0 χ 0 → γγ (or γZ) at (or very close to) E γ = m χ , where m χ is the mass of the dark matter particle, to be determined.
In particular, the paradigmatic WIMP with mass in the 100 GeV to 10 TeV range and electroweak charge is expected to be observed or ruled out by the Cherenkov Telescope Array (CTA) [1] under construction even under conservative assumptions on astrophysical uncertainties, especially due to the dark matter density profile near the Galactic center. Precise theoretical computations of the photon yield from DM annihilation are therefore well motivated.
Recent theoretical work has focused on two aspects of the problem. First, for dark matter annihilation into energetic particles, electroweak Sudakov (double) logarithms O((α 2 ln 2 (m χ /m W )) n ) are large and should be summed to all orders [2][3][4][5], in addition to the summation of ladder diagrams known as the Sommerfeld effect. Second, since γ-ray telescopes do not measure two photons from a single annihilation in coincidence, the observable is not χ 0 χ 0 → γγ (or γZ) but rather the semi-inclusive single-photon energy spectrum γ + X, where X denotes the unidentified other final state particles. Although the leading term in the perturbative expansion of the semi-inclusive annihilation rate arises from the two-body final states γγ, γZ, the logarithmically enhanced terms differ in higher orders and this affects their resummation [6][7][8]. It has been shown, both for the exclusive γγ annihilation rate [5], as well as for the semi-inclusive rate at narrow energy resolution (as defined below) [7], that resummation with NLL' accuracy, which combines the full one-loop calculations with next-to-leading logarithmic resummation provides precise results for the photon rate with uncertainties around 1%.
The resummation of the semi-inclusive spectrum is performed for the primary photon energy spectrum d(σv rel )/dE γ of the DM pair annihilation cross section multiplied by the relative velocity of the annihilating particles. While in forecasts for the rate observed by a specific telescope, the spectrum will have to be smeared with an instrument-specific resolution function of some width E γ res in energy, the expected impact and accuracy of the theoretical prediction can be equally discussed for the spectrum integrated over the energy interval E γ res from its kinematic endpoint: The endpoint-integrated spectrum depends on the three scales m χ , m W (representative of electroweak scale masses), and E γ res . We consider TeV scale dark matter, hence the hierarchy m W m χ is always assumed. The details of the resummation of electroweak Sudakov logarithms near the endpoint, E γ res m χ , differ according to the scaling of E γ res and m W with respect to each other. We distinguish the following three regimes: The wide resolution regime was considered in [6,8] and resummed at the NLL order. Due to the double hierarchy m χ E γ res m W a two-step procedure applies to simultaneously sum the unrelated large logarithms of m χ /m W and E γ res /m W . This procedure requires large dark matter masses to satisfy both hierarchies. Resummation of electroweak Sudakov logarithms for the narrow resolution case was accomplished in [7] at the NLL' order. The intermediate resolution regime has not been considered up to now.
In the present paper we close this theoretical gap. We develop the effective field theory (EFT) for the intermediate resolution regime and sum the electroweak logarithms at the NLL' order. We show that the result can be smoothly joined to the narrow resolution regime to provide a precise prediction of the photon energy spectrum near m χ in the entire region from the line signal (E γ res = 0) to E γ res ≈ 4m W . We also provide details and derivations for the narrow resolution regime not given in the letter [7].
The intermediate resolution regime is relevant to present and upcoming DM searches. For example, assuming the regime to apply to E γ res in [m W /4, 4m W ] the energy resolution of the H.E.S.S. experiment E γ res /E γ ≈ 10% [10] implies that dark matter masses in the range 200 GeV to 3.2 TeV are covered by the intermediate resolution calculation. For the CTA experiment, we obtain the power-law fit E γ res /E γ = 0.0915 (E γ /TeV) −0.347 from Figure 11 of [9] in the range of photon energies of interest, which is shown as the dash-dotted line in Figure 1 together with the unapproximated resolution (solid line). The horizontal band (dark-grey/red) represents the region of applicability of the intermediate resolution regime, which extends to 6.8 TeV for the CTA experiment. Thus, the intermediate resolution calculation covers the mass values where the thermal relic density of the pure Higgsino (electroweak doublet) and pure wino (triplet) models agrees with the observed relic density.
The outline of the paper is as follows. In Section 2 we discuss the momentum modes and effective Lagrangians relevant to the problem and derive a factorization formula for the photon energy spectrum valid when m χ − E γ = O(m W ), which corresponds to the intermediate resolution regime. We also discuss the modifications that apply to narrow resolution [7]. In Section 3 we calculate the hard, jet and soft functions that appear in the factorization formula at the one-loop order, as required for NLL' resummation, provide the renormalization group equations (RGEs), and solve them with corresponding accuracy. These calculations are performed for the specific case of the pure wino model, which corresponds to the Standard Model extended by an SU(2) triplet with zero hypercharge, of which the electrically neutral member is the dark matter particle. In the subsequent Section 4 we show our main result, the resummed endpoint-integrated photon spectrum in the range of dark matter masses of interest and for various E γ res . We match the intermediate resolution calculation to the narrow resolution one from [7] and find very good agreement. We pursue and explain this numerical observation in Section 5 by expanding analytically the resummed expressions to the two-loop order and comparing the logarithmic and constant terms. We conclude in Section 6. A series of appendices collects additional technical details on soft and jet function integrals, including the narrow resolution case and the treatment of the Z-boson resonance, the RGE invariance check for the narrow resolution case, and the complete analytic expressions for the expansion of the resummation formula to the two-loop order.

Factorization of the energy spectrum
The annihilation cross section of TeV scale dark matter with electroweak charges can be strongly modified by the Sommerfeld effect [11] due to non-relativistic scattering of the dark matter particles before they annihilate. This effect is well understood. Our concern are the large electroweak logarithms in the annihilation rate χ 0 χ 0 → γ +X, when the photon energy is close to maximal, E γ ∼ m χ m W , more precisely m χ − E γ ≤ E γ res m χ . The observation of a photon with this energy implies that the unobserved final state X is "jet-like" with small invariant mass m X = 4m χ E γ res . The logarithmic enhancements of such final states are caused by soft and collinear physics relative to the large scale m χ . In a systematic expansion in m W /m χ , where E γ res is parametrically related to m W , m χ as above, the χ 0 χ 0 → γ + X process is separated into a hard annihilation process and the low-scale initial-and final-state dynamics, which is described by suitable effective Lagrangians valid at scales µ m χ . We assume that the DM particle is the component of an SU(2) multiplet of the electroweak interaction, which remains electrically neutral after electroweak symmetry breaking. Since m χ m W , this is always a good approximation at leading order in an expansion in m W /m χ unless there are two nearly degenerate heavy multiplets, such that large mixing can occur. For definiteness, we assume (as in [7]) that χ a , a = 1, . . . , 2j + 1, is a 2j + 1 dimensional isospin-j SU(2) multiplet of Majorana fermions with integer j (thus hypercharge vanishes). The essence of the derivation of the factorization formula below does not rely on these assumptions.

Effective Lagrangians and annihilation operators
After integrating out virtualities of order m 2 χ , the short-distance part of the annihilation process is represented by an operator that destroys the two DM particles at a single point, and a set of collinear and anti-collinear fields along opposite light-like directions starting from this point, which describe the energetic particles in X and those that convert to the observed photon. We refer to the direction n µ − of the jet X as "collinear". The direction of the photon momentum defines the "anti-collinear" direction, p µ γ = E γ n µ + . The reference vectors satisfy n 2 − = n 2 + = 0, n + · n − = 2. A general momentum is written in components as k µ = (n + k, n − k, k ⊥ ), such that for collinear momenta n + · k n − · k, vice-versa for anti-collinear momenta.
The low-energy dynamics of the intermediate resolution case is described by nonrelativistic effective field theory [12] and soft-collinear effective field theory (SCET) [13][14][15]. The kinematics of the annihilation process considered here is a mixture of an inclusive process in the collinear direction, also called a SCET I problem, and an exclusive final state of the SCET II type in the other direction with the added complication of electroweak symmetry breaking and gauge boson masses. The effective Lagrangian must describe the interactions of the relevant modes with momentum scaling hard-collinear (hc) : k µ ∼ m χ (1, λ, √ λ) collinear (c) : k µ ∼ m χ (1, λ 2 , λ) anti-collinear (c) : k µ ∼ m χ (λ 2 , 1, λ) soft (s) : k µ ∼ m χ (λ, λ, λ) Hard modes with momentum k µ ∼ m χ (1, 1, 1) are integrated out into matching coefficients and are no longer part of the effective Lagrangian by construction. The power counting parameter is the small ratio λ = m W /m χ . Compared to the narrow resolution case [7], an additional hard-collinear mode is needed to describe the unobserved final state X with hard-collinear virtuality of order m 2 χ λ = m χ m W . On the other hand, the effective theory for the wide resolution case [6,8] requires a yet more numerous set of modes to account for the independent scales E γ res and m W . This set collapses to the one above when E γ res is set parametrically to m W . The leading hard annihilation processes are those into two energetic final-state particles. Adding another collinear or anti-collinear field to the primary annihilation vertex, implies a suppression by at least one power in λ due to the scaling of the fields in the effective Lagrangian. In this work as well as in all previous works on electroweak Sudakov resummation for dark matter annihilation, the aim is to sum logarithms of m W /m χ at leading power in the expansion in λ. Power-suppressed effects in λ = m W /m χ are systematically neglected in this treatment.
A general consequence of neglecting power corrections is that (anti-) collinear fields must preserve their identity while emitting soft radiation. Since the energetic final state in the anti-collinear direction consists of a single photon with nearly maximal energy, which hence cannot be generated from an energetic fermion or Higgs boson, the leadingpower operators for the hard annihilation process contain a single anti-collinear SU (2) gauge field. The collinear part of the operator must then also be an SU(2) gauge field by gauge invariance, because the non-relativistic initial state consists of a dark-matter two-particle state with vanishing hypercharge and colour, and integer weak isospin. It is therefore not possible to combine the anti-collinear SU(2) gauge field with any other single Standard Model (SM) field to form these quantum numbers, except with another SU(2) gauge field. The hard annihilation process is reproduced by the effective Lagrangian with operators of the form Here χ v is a two-component non-relativistic spinor field in the SU(2) weak isospin-j representation, χ c v = −iσ 2 χ * v = − χ * v (with the antisymmetric 2 x 2 matrix with 12 = 1) the charge-conjugated field, and A A ⊥c,µ (A B ⊥c,ν ) the collinear (anti-collinear) SU(2) gauge field of soft-collinear effective theory. The definitions will be given below. Fields without position arguments are evaluated at x = 0. The operator is non-local, since (anti-) collinear field operators are integrated along the light-cone of the respective direction with the coefficient functionĈ i . 1 The spin matrix Γ µν i is contracted with the two-spinor indices of the DM fields (not written explicitly) and the Lorentz index of the gauge fields. Similarly, the SU(2) group index matrix T AB i is contracted with the two isospin-j representation indices of the DM fields (not written explicitly) and the adjoint index of the gauge fields. The operator basis is given by the list of distinct spin-and group-matrix structures.

Non-relativistic dynamics
For energies below m χ but above m W the DM interactions are described by the standard non-relativistic Lagrangian with D µ = ∂ µ − ig 2 A C µ T C the SU(2) covariant derivative, T C , C = 1, 2, 3, the SU(2) generators in the isospin-j representation of the DM field, and g 2 the SU(2) gauge coupling. The Lagrangian can be extended to include interactions suppressed by powers of p/m χ . Since the largest non-relativistic momentum scale is m W , they correspond to power corrections, which are neglected.
The non-relativistic Lagrangian describes the soft, potential and ultrasoft modes (see, for example, [16]) of the non-relativistic DM and light SM fields. The soft modes can be integrated out from the non-relativistic Lagrangian in straightforward analogy with heavy quark anti-quark systems in non-relativistic QCD. Together with the potential modes of the light particles, they generate instantaneous but spatially non-local interactions between the DM fields, that is, DM potentials. The effective Lagrangian for the remaining potential modes of the DM field and the ultrasoft modes of the light fields, is the potential-non-relativistic dark matter Lagrangian [17]. At leading power, We indicated explicitly the space-time arguments of the fields to highlight the nonlocality of the potential interaction and the fact that the ultrasoft gauge field in the covariant derivative D 0 is multipole-expanded around x = 0. Note that the covariant derivative is now the one with respect to the unbroken electromagnetic gauge symmetry, since ultrasoft light fields with momentum k ∼ m χ λ 2 ∼ m 2 W /m χ exist only for fields with masses much smaller than m W . The electroweak gauge bosons no longer appear as dynamical fields in PNRDM effective theory.
The soft modes of the light particles have virtuality of order m 2 W . Hence, in matching NRDM EFT to PNRDM EFT, the masses of the electroweak gauge bosons and of the top quark and Higgs boson cannot be neglected. The potential V {ij}{kl} (r) depends on these masses, resulting in Yukawa (electroweak gauge bosons, Higgs bosons) and Coulomb potentials (photons). Furthermore, the components of the original isospin-j DM multiplet acquire slightly different masses after electroweak symmetry breaking. The Lagrangian above uses δm i = m i − m χ 0 ≥ 0, where m i is the mass of eigenstates labelled by i. Since L PNRDM is no longer invariant under the SU(2) L ×U(1) Y gauge symmetry and calculations are carried out in broken theory, we express it in terms of mass eigenstate fields χ vi rather than the fields χ va of the SU(2) multiplet.
While soft subgraphs not connected to the annihilation vertex generate potential interactions, soft momentum running through the annihilation vertex, dresses the nonrelativistic fields in the operators (5). Since the leading soft interaction is of the eikonal type, this dressing takes the form of a Wilson line. This is seen most easily by noting that the temporal soft gauge-field coupling in the covariant derivative D 0 in (6) can be removed by the field redefinition where the soft Wilson line Y v (x) is defined as the path-ordered exponential with T C the SU(2) generators in the spin-j representation and v µ = (1, 0). 2 Dropping the superscript (0), the non-relativistic part of the operators (5) takes the form after the field redefinition.
The main use of the PNRDM Lagrangian after the field redefinition is related to the computation of the Sommerfeld effect. In this context it is convenient to write the sums over the field indices in terms of a sum over the composite indices I = {ij} and K = {kl} of two-particle states according to the bound-and scattering-states of the Schrödinger problem for the relative coordinate. For example, for the triplet wino model, the index i takes the values 0, +, − corresponding to the electric charge of the DM mass eigenstates. The index I runs over the nine two-particle values 00, +−, −+, 0+, +0, 0−, −0, ++, −−, ordered by the modulus of the electric charge. Since electric charge is conserved, it is sufficient for the computation of the χ 0 χ 0 annihilation rate to solve the Sommerfeld problem in the charge-0 sector of the two-particle states. For the SU(2) j = 1 triplet the fields are related to the mass basis by χ ± = (χ 1 ∓ iχ 2 )/ √ 2, χ 0 = χ 3 . The two-particle states are related by where the 9 × 9 matrix K ab,I can be read off from  While the specific form of the K-matrix depends on the SU(2) representation of the DM field, the formalism is general. 2 The field redefinition is analogous to but not the same as the field redefinition discussed in [18]. There ultrasoft gauge bosons are decoupled from potential fields in potential non-relativistic EFT. Here the decoupling refers to soft gluons in NRDM EFT. In both cases the soft Wilson lines can and must be evaluated at the multipole-expanded position x 0 = (t, 0), since the three-momentum of the gauge boson does not enter the virtual heavy particle propagator in the leading-power approximation.
From this point on the non-relativistic part of the problem follows the discussion of the computation of the Sommerfeld effect for an arbitrary set of heavy fermions nearly degenerate with the DM particle, developed for the general minimal supersymmetric SM in [17,[19][20][21]. The framework described in these papers in turn generalizes the original work [11,22] to mixed DM states and reformulates it in the DM EFT context. For the case of the pure wino triplet, the result is identical to the original treatment, but can in principle be extended to systematically include radiative and velocity corrections to the Sommerfeld effect. In this paper, however, as in all previous studies, the Sommerfeld effect is computed only at leading order in the PNRDM EFT.
What will appear in the factorization formula is the matrix element of non-relativistic annihilation operators of the form which arise from (10) after squaring the amplitude. The matrix Γ = 1 or σ i (Pauli matrix) operates on the spinor index of χ v , depending on whether the fermion bilinear destroys or creates a spin-0 or spin-1 state. By assuming that the spin matrices in the two bilinears are the same, we implicitly make use of the fact that the potential V {ij}{kl} (r), while being spin-dependent, does not change the spin of the incoming two-particle state before it annihilates. The NRDM EFT matrix element of the above operator in an incoming χ i χ j DM state with some relative velocity v rel , orbital quantum number L = 0 (S-wave) and total spin S is given by (no sum over i, j) [20] where ψ (L,S) e 1 e 2 , ij is the χ e 1 χ e 2 -component of the scattering wave function for the incoming χ i χ j state, evaluated for zero relative distance and normalized to the free scattering solution, that is ψ (L,S) eae b , ij → δ eai δ e b j in the absence of interactions. The symbols ξ i , ξ j in the second line of (14) denote the Pauli spinor of the incoming particles χ i and χ j , and . . . stands for the trace in spin space (spin sum). The multi-component wave function ψ e 2 e i ,ij accounts for the potential interactions of the incoming χ i χ j state, which couple it to all possible intermediate two-body states e 2 e 1 with the same charge, spin and orbital angular momentum. Both wave-function components e 1 e 2 and e 2 e 1 contribute to the matrix-element of the operator χ c † e 2 χ e 1 . For an operator with quantum numbers L and S, there is a relative sign (−1) L+S between the two components. The above-defined ψ (0,S) e 1 e 2 , ij is related via to the coordinate-space scattering wave-function [ψ E (r )] I,ij at the origin, which in turn can be obtained directly from the matrix-Schrödinger equation with the potential (7), now including the mass splitting between the mass M I of the twoparticle state I and the mass of the χ 0 χ 0 state via V IK → V IK + δ IK (M I − 2m χ ). The energy E is fixed through the relative velocity of the initial state, and the label ij refers to the fact that this equation should be solved with the initial condition corresponding to the particular incoming two-particle state ij. µ I is the reduced mass of the two-particle state I, which can be set to m χ /2 in the leading-order treatment of the Sommerfeld effect. We refer to [20] for further details and the methods employed to solve this equation.
For the task at hand, we focus on the initial state I = ij = 00. The other two-particle states appear only as virtual states in the ladder diagrams summed by the Schrödinger equation. 3 Due to electric charge conservation, the potential is block-diagonal and it is sufficient to solve (16) in the charge-0 sector, which, for the wino example, consists of I = 00, +−, −+. The description in terms of three two-particle states is convenient, since the framework can be formulated without additional rules for the construction of the potential for different S and L. The (anti-) symmetrization is encoded automatically in the (anti-) symmetry of the operator and its short-distance coefficient. However, the description is redundant, since the fermion bilinear with +− fields is identical up to a possible sign to the one with −+. It is customary in the discussion of the Sommerfeld effect to reduce the basis of two-particle states to non-identical ones (six instead of nine, for the triplet model, and two instead of three for the charge-0 sector). In the following we adopt this convention. Specifically, for the wino (triplet) model, I shall then refer to 00, +− only. The precise relation between the two formulations, referred to as method-1 and method-2, respectively, is explained in [20], including the explicit forms of the potential and tree-level short-distance annihilation coefficients in both methods. Irrespective of the method, the Sommerfeld factors are defined as The discussion up to now ignored the coupling of ultrasoft photons to the charged members of the DM multiplet through the electromagnetic covariant derivative in (7). This is justified, since the field redefinition mentioned in footnote 2 removes the ultrasoft photon field from the Lagrangian at the expense of modifying the DM fermion bilinear as where S vi is an electromagnetic time-like Wilson line corresponding to the charge of the field χ vi in I = {ij}. In the charge-0 sector, the charges of the fields i, j add to zero, which implies S vi S vj = S vi S † vi = 1, reflecting the well-known fact that photons with wave-length much larger than the size of the system only couple to the total charge (here, zero) of the system.
Finally we note that the factorization of non-relativistic dynamics from the soft and collinear dynamics of the final state is independent of the photon energy resolution as defined above, at least to the accuracy considered here. The key requirement is the decoupling of soft and ultrasoft interactions from the ladder diagrams that build up the Sommerfeld effect, which holds because soft gauge bosons throw potential DM propagators off-shell, and because ultrasoft photons do not interact with the electrically neutral two-particle state.

Soft-collinear dynamics
The characteristics of the final state is an energetic photon, whose momentum is balanced by a jet of unobserved particles. The low-energy physics of energetic objects with small invariant mass interacting with soft modes is described by SCET. Its application to electroweak Sudakov situations was first discussed in [23,24] for the production of two particles with electroweak charges in a high-energy collision. As in the non-relativistic sector, one needs two different effective Lagrangians depending on whether the virtuality is much larger than or of order m 2 W . Although in higher orders all SM fields are present in the collinear and soft interactions, we restrict ourselves to the gauge boson Lagrangian, since the gauge boson SCET fields appear directly in the annihilation operators, and since the discussion of fermions is fairly standard from QCD applications of SCET. The SCET Lagrangian consists of where L soft is the purely soft field Lagrangian that takes the same form as the corresponding SM Lagrangian except that all fields are assumed to be soft. The collinear Lagrangian at leading power is where g 2 F µν c = i [D µ , D ν ] as usual, but the collinear SU(2) covariant derivative is given by We included the collinear Higgs doublet field ϕ c for later purposes. At leading power, soft-collinear interactions involve only the n − A s projection of the soft gauge field. Moreover, the soft gauge field is evaluated at the multipole-expanded position x µ − + x µ ⊥ with x µ − = (n + x) n µ − /2, reflecting the fact that the n + k component of the soft momentum can be neglected relative to the corresponding large component of hard-collinear and collinear momentum. The Lagrangian above accounts for collinear modes of both, the hard-collinear and collinear type and is formulated in the unbroken phase of SU(2) gauge symmetry, since this is relevant to the hard-collinear fields of virtuality m χ m W . The anti-collinear Lagrangian L c is the same up to the interchange of n + ↔ n − . The above expressions should be amended in an obvious manner to include the gauge field for the U(1) hypercharge interaction and its coupling to the Higgs field.
The SCET Lagrangian enjoys separate collinear, anti-collinear and soft gauge symmetries [25]. It is convenient to express the collinear Lagrangian in terms of manifestly collinear-gauge invariant collinear fields 4 where the collinear Wilson line is defined as In terms of these, together with the collinear Lagrangian is expressed as In this form it is apparent that the collinear gauge field degrees of freedom are represented by the transverse fields, since n + A c = 0 from (22), while n − A c can be eliminated using the gauge-field equation of motion (see, for instance [26], Appendix B for the operator equation in QCD). At scales µ m χ there are no interactions between collinear modes of different directions, as well as between collinear and non-relativistic DM modes, since this would result in hard virtualities, which are already integrated out into the short distance coefficients of annihilation operators. However, they all interact with each other through the soft gauge fields. As was the case for the non-relativistic DM field, the soft gauge field can be decoupled from the hard-(anti-)collinear fields through the field redefinition 5 with [14] Y 4 The following construction can be extended to the hypercharge gauge field by introducing an abelian U(1) collinear Wilson line. Since the Higgs field carries hypercharge, the collinear-invariant field includes the hypercharge Wilson line. 5 In the coupling to hard-collinear (as opposed to collinear) fields, the argument of the soft field in the covariant derivative (21) can be set to x − at leading-power accuracy, since the transverse momentum of the soft mode is negligible compared to the hard-collinear one. This simplification is required, so that the field redefinition in the following equation removes the soft field from the covariant derivative on hard-collinear fields. Also note the order of adjoint indices in (26). With the alternative definition Here the SU(2) generator T D refers to the adjoint representation, (T D ) BC = −i DBC , in case of the adjoint gauge field, and the fundamental representation for the analogous decoupling transformation of the (anti-) collinear Higgs fields. As a result of this field redefinition, soft Wilson lines appear in the annihilation operator.
In the intermediate resolution case, the virtuality of the unobserved jet is not resolved by the measurement below the hard-collinear scale m χ m W . The dynamics of this jet is described by hard-collinear modes. On the other hand the scale for the anti-collinear direction of the observed photon is set by the virtuality m 2 W of the collinear electroweak gauge bosons, whose masses cannot be neglected. The photon "jet-function" as well as possible mass effects within the hard-collinear jet must be computed with the SCET Lagrangian for the (anti-) collinear modes of the massive electroweak gauge bosons and the photon after electroweak symmetry breaking. The gauge boson mass term follows from the Higgs covariant kinetic term in (25) from which shows explicitly that the mass term arises only for the transverse field. The collinear gauge field Lagrangian for virtualities of order m 2 W , ignoring now the Higgs field, reads The hypercharge interaction should be added in the standard way. It is then convenient to express the Lagrangian in terms of mass eigenstates fields. The collinear fields no longer interact with soft fields, as discussed above, but the interaction with ultrasoft fields is still present. However, only fields with masses much smaller than m W can be ultrasoft, and the leading-power ultrasoft interactions are included through electromagnetic covariant derivatives acting on the electrically charged electroweak gauge fields with covariant derivatives defined as in (21), except that now n − A s (x − ) refers the ultrasoft photon field only.

Annihilation operator basis
The annihilation operators (5) and their short-distance coefficients do not depend on the photon energy resolution as long as E γ res m χ . Two relevant operators have been identified in previous work [4,5,7], but the arguments that these two operators form a complete basis to all orders in perturbation theory have not been explicitly provided. In the following we use symmetries to reduce the possible operators to with the spin matrix in d space-time dimensions given by (Conventions v µ = (1, 0, 0, 0), n µ ± = (1, 0, 0, ∓1), n = (0, 0, 1), m, n = 1, 2, 3, 0123 = −1 are used.) We then show that the third operator does not contribute when the detected gauge boson is a photon.
We have already shown that the collinear and anti-collinear field must each consist of a single SU(2) gauge field. We therefore start from the general form (5) and note that the two DM fields must couple to an operator with SU(2) isospin 0, 1 or 2. Thus, the group-index matrix must be from where T A are the SU(2) generators in the isospin-j representation.
Turning to the spinor and Lorentz indices, the two spin-1/2 DM fields can couple to spin-0 or spin-1. In the first case the implicit pair of two-spinor indices of Γ µν must be of the form δ αβ . The spin-1 structure is the vector of Pauli matrices (0, σ) or [σ ρ −(v·σ)v ρ ] αβ . For the spin-0 case, noting that µ, ν are transverse indices, we obtain two different Γ µν by multiplying with defined in (33). For spin-1, the three independent combinations can be formed. Here v ρ (σ ρ − (v · σ)v ρ ) = 0 was used to reduce a number of further structures to the given ones. Together with the three independent SU(2) structures, this results in six spin-0 and nine spin-1 operators. The final state of two gauge bosons must respect Bose symmetry, hence the operator has to be symmetric under the simultaneous exchange of all labels, c ↔c, n + ↔ n − , A ↔ B, µ ↔ ν. The admissible structures are therefore the product of T AB 1 , T AB 2 (T AB 3 ) with the symmetric (antisymmetric) tensors from (36), (37), which leaves three spin-0 and four spin-1 operators.
The DM gauge interaction conserves CP symmetry and consequently parity for Majorana fermions. Since χ c † v χ v (χ c † v σχ v ) has negative (positive) parity, this excludes g µν ⊥ in (36) and all except the first structure in (37), resulting in the two spin-0 operators O 1,2 and one spin-1 operator O 3 as given above.
In the process χ 0 χ 0 → γ + X, the assumption of a single photon in the anti-collinear final state implies that the SU(2) index B in O i must necessarily be B = 3. For the operator O 3 , the index C in χ c † v T C χ v must then be C = 1 or 2, in which case the bilinear cannot annihilate an electrically neutral two-particle state. Thus, O 3 does not contribute to the annihilation into γ + X when the photon is required to have nearly maximal energy. We note that the remaining two operators are both spin-singlet, so the dominant short-distance annihilation process occurs in the 1 S 0 configuration.
When the matching calculations are done with dimensional regularization, the question arises whether the above operator basis is complete in d dimensions or whether it has to be complemented by evanescent operators, which vanish in d = 4. We find that no evanescent operators arise. The d-dimensional basis is given by the same O 1,2 except that the first form on the right-hand side of (33) should be used for the spin matrix rather than the four-dimensional expression µν ⊥ . To see this we note that an arbitrary full theory diagram in the calculation of the hard matching coefficients contains a single string of Dirac matrices of the form 6 with indices µ i that can be contracted among each other, with the vectors v or n ± and with the two transverse polarization vectors ε c⊥ , εc ⊥ of the external gauge boson lines. To obtain an S-wave annihilation operator in the non-relativistic EFT, N must be odd. By systematically exploiting the on-shell condition / vu(p) = u(p), the relations n − = 2v − n + and n + · ε c⊥ = n + · εc ⊥ = 0, which imply {/ n + , / ε c⊥ } = {/ n + , / εc ⊥ } = 0 the string can be reduced to the two structures After expressing the Dirac spinors in terms of non-relativistic two-spinors, the first structure corresponds to the spin matrix of O 3 , and the second to Γ µν . The operator basis holds for any integer isospin-j DM multiplet with vanishing hypercharge. The coefficient functions C 1,2 of O 1,2 and their renormalization group evolution (RGE) to scales µ m χ can be found in [7] and we refer to this paper and Section 3.1 below for the detailed expressions.

Factorization
In this section we derive the factorization formula for the photon energy spectrum for intermediate photon resolution. We also comment on the modifications for the narrowresolution case, in this way providing the derivation of this case omitted in [7]. To guide the reader let us preview here the main result of this section by giving the equation that will subsequently be used for the calculation of the resummed spectrum at the NLL' order.
Independent of the resolution we can represent the energy spectrum in the form where the sums over I, J run over all electrically neutral two-particle states that can be formed from the 2j + 1 single-particle states of the electroweak DM multiplet, and the sums over i, j refer to the two operators O 1,2 . The expression after the first equality expresses the factorization of the Sommerfeld enhancement factor from the remainder of the process. S IJ is the same Sommerfeld factor as usual, except that the tree-level short-distance annihilation matrices are replaced by matrices Γ IJ (E γ ), which include electroweak Sudakov resummation and other radiative corrections up to the specified accuracy.
The various functions will be defined below. Some of them require rapidity regularization in addition to the conventional dimensional regularization, resulting in a dependence of the renormalized function on the rapidity factorization scale ν in addition to the dimensional regularization scale µ. We note the convolution of the jet function J int for the unobserved final state X with a soft function W , which accounts for radiation of soft electroweak gauge bosons and other soft particles into the final state, and virtual corrections. We can compare this to the corresponding formula for narrow resolution, The main difference is that the smaller invariant mass of the final state X forbids soft real radiation. The soft effects are purely virtual, and appear at the amplitude level in the factors D.
The starting point for the derivation is the general expression for the initial-state spin-averaged and final-state spin-summed annihilation cross section The sum-integral symbol implies a sum over all kinematically allowed final states X with total momentum p X and the phase-space integral over the final-state momenta. Summation over spins is understood for the initial and final state and the overall factor 1/4 accounts for the initial-state spin average. In the center-of-mass frame the initialstate momentum is p χχ = (2m χ + E)v. T i→f is the T-matrix element for the transition. After integrating out the hard momentum modes, the T-matrix element is non-vanishing only if it involves the effective interaction (4), and we can write We have split the sum over X into a sum over (hard-) collinear particles X c and soft particles X s . The matrix element is to be evaluated in non-relativistic and soft-collinear EFT. The factor 2m χ arises from the non-relativistic normalization of external DM state. After the field redefinitions (8), (26) that decouple soft gauge bosons from the (hard-) collinear, (hard-) anti-collinear and non-relativistic fields the operators are We now use the symbol Y ± (x) to denote SU(2) Wilson lines in the adjoint representation and recall that fields and Wilson lines without space-time argument are evaluated at x = 0. Because the different types of fields no longer interact, we can factorize the matrix element into Translation invariance implies where p Xc is the total four-momentum of the collinear final state, which allows us to perform the s, t integrations in (44) and express them in terms of the momentum-space coefficient function Up to power-suppressed corrections n − p γ = 2E γ ≈ 2m χ , n + · p X ≈ 2m χ . We therefore define and these are given in [7] and below in the one-loop approximation required for NLL' accuracy.
We write the four-momentum conservation delta-function in (43) as the space-time integral of the exponential, insert the factorized matrix element γ(p γ )X c X s | O i |[χχ] 00 into (44), and the square of the resulting expression for the T-matrix element in (43). In this way, we obtain We use translation invariance again to absorb the exponentials e −ip Xc ·x , e −ip Xs ·x into a shift of position 0 to x in the first halves of the collinear and soft matrix elements, after which the sums over complete sets of collinear and soft intermediate states can be done. We also use (14), (17) to express the non-relativistic matrix element in the form with ξ 0 the spinor of an external χ 0 field (with the two orientations ↑, ↓). The Sommerfeld factor is a function of the small kinetic energy E of the DM two-particle state. For annihilation in the present Universe E is much smaller than any other energy scale in the problem. After factoring the Sommerfeld effect, E can be neglected in the other parts of the calculation, that is, we set p χχ = 2m χ . After these manipulations in (47), by comparing to (40) we can read off the quantity γ ij IJ (E γ ): introducing the soft operator The last three factors in the above equation define, in order, an anti-collinear, hardcollinear and soft function, as follows.

Definitions for the intermediate resolution case
Photon collinear function The "jet" function for the exclusive anti-collinear photon state is defined by the squared matrix element We have made the sum over photon polarizations explicit. Obviously, only Z 33 γ is nonvanishing. From (22) and (23) is follows that Z 33 γ /ŝ 2 W can be interpreted as the on-shell photon field renormalization constant in anti-collinear light-cone gauge n − · Ac = 0. Z 33 γ depends on the electroweak scale masses m W , m Z , m H and m t of the SM particles, the dimensional regularization scale µ and a rapidity regulator scale ν, since the factorization formula involves the separation of regions (here anti-collinear and soft) with equal virtuality but parametrically different n ± momentum components.
Unobserved-jet collinear function The jet function pertaining to the inclusive (unobserved) collinear final state is defined as The jet function is defined in SCET I in terms of the hard-collinear gauge field. It depends on the hard-collinear scale through the invariant mass squared p 2 of the final state X, but also on the scale m W through the electroweak scale masses of the particles inside the jet. The jet function as defined above is therefore still a two-scale object, which can be further factorized into a hard-collinear and collinear function [6]. Up to power corrections of order m 2 The hard-collinear matching coefficient J int (p 2 ) can be computed in the theory with unbroken electroweak gauge symmetry in close analogy with the standard gluon jet function in QCD. It depends on the renormalization scale µ, but does not require rapidity regularization. The collinear mass-jet function J XV m (m W ) is momentum-independent, but can depend on both µ and the rapidity regulator. However, we find that at tree-level and at the one-loop order, the collinear mass function is trivial, that is It is plausible that this result holds to any order in the coupling, since the observable is not sensitive to the internal jet structure. 7 We shall make use of this simplification in deriving (41).

Soft function
The sum over the soft final state in (52) is the unit operator, which allows us to define the soft function in momentum space, We also define the integrated soft function where ω = n − · k, and then the SU(2) index-contracted soft function The soft functions must be calculated in the broken SU(2) theory and depend on the electroweak masses of the SM particles. They also depend on the renormalization scale µ and the rapidity regularization scale ν.

Derivation of the final formula
With the above definitions of the collinear and soft function, we can rewrite the corresponding terms in (52), where in passing from the second to the third line we used There is no dependence on the direction of the photon momentum, hence we can perform the photon phase-space integral in (52), to obtain This equation represents the factorization formula for the intermediate resolution case in its general form. To obtain (41), we note that both operators involve the same spin matrix (33), that is Γ µν We then use the property (57), which allows us to replace Finally, we switch from method-1 to method-2 (see the discussion before (17)) and sum only over distinguishable two-particle states I, J. As discussed in [20], this implies certain replacement rules for the potential used in the computation of the Sommerfeld effect and the annihilation matrix Γ IJ , which introduces the factor 1/( √ 2) n id in (41), where n id = 0, 1, 2 depending on how often the two-particle state 00 appears in the index pair IJ. The objects in the factorization formula are assumed to be evolved from their natural scales, where they exhibit no large logarithms, to a common scale in µ and in ν with the renormalization group equations discussed in the following section. The evolution factors accomplish the desired resummation of large logarithms.
Let us comment on the treatment of ultrasoft modes that we did not mention in the derivation. After the decoupling of soft modes from the (anti-) collinear and nonrelativistic fields, all of them still interact with ultrasoft modes. In writing the various sectors in a factorized form, we implicitly made use of the fact that an ultrasoft Wilson line field redefinition decouples ultrasoft interactions from (anti-) collinear modes at leading power. This introduces multiple convolutions from the different sectors with an ultrasoft function. We omitted this ultrasoft function in the above discussion, since it is actually absent due to the electric charge neutrality of the initial state and the anticollinear photon final state. The soft and the hard-collinear final state, however, are not necessarily electrically neutral. However, all momentum components of an ultrasoft mode are small compared to the corresponding momentum component of a soft or hardcollinear mode, such that in leading power, the ultrasoft momentum transfer to the soft or hard-collinear function can be neglected. This eliminates the possibility of a non-trivial convolution and allows us to ignore the ultrasoft mode.

Modifications for the narrow resolution case
Following the above line of argument, we derive the factorization formula for the narrow resolution case stated in [7] and written in (42) in present notation. There is no change to the discussion of the non-relativistic and photon jet function, but for the unobserved jet function and soft function, the following differences need to be noted.
The narrow resolution jet function has the same definition as (55), except that now the gauge field is collinear rather than hard-collinear. Consequently, there is no further factorization. Since there is no soft radiation into the final state, the collinear function must be charge-neutral which selects the X = V = 3 component of the jet function. The narrow resolution jet function depends on the collinear invariant mass squared p 2 ∼ m 2 W and the electroweak scale particle masses of the SM. It further depends on the renormalization scale µ and, contrary to the intermediate resolution hard-collinear jet function, also on the rapidity scale ν. The different structure of rapidity logarithms for the two cases is matched by different rapidity logarithms in the soft function. Details on the calculation of the narrow resolution jet function are provided in Appendix B.2.
The small energy resolution E γ res ∼ m 2 W /m χ forbids soft radiation into the final state, hence the soft factors are defined at the amplitude level, rather than for the square of the amplitude as above. Technically, the sum over the soft final state in (52) is empty, such that where D i I,V W is defined as the vacuum matrix element of the soft operator (53). The photon jet function selects Y = W = 3, and the unobserved jet function selects X = V = 3, which implies that only the single SU(2) component D i I, 33 of the soft amplitude is needed. As in the intermediate resolution case the soft function must be computed in the broken theory. The inclusive nature of the soft function for intermediate resolution entails a partial cancellation of infrared singularities between virtual and real contributions, which does not happen for the narrow resolution case, where real contributions are absent. This also changes the structure of the rapidity evolution factor, since the narrow resolution soft function couples to the rapidity evolution of the collinear and anti-collinear sector, while only the latter has rapidity divergences for intermediate resolution. This explains the differences between the two factorization formulas (42) and (41). 8 ization theorems for the semi-inclusive photon spectrum in DM annihilation derived in the previous section. Furthermore we will show the consistency of the renormalization group and discuss different resummation schemes.
The hard functions have been computed for an electroweak DM with any integer isospin j. The (anti-) collinear functions for the photon and for the unobserved jet triggered by an electroweak gauge boson are universal. The soft function given below is specific to the triplet (wino, j = 1) DM model, which is the focus of this work.

Hard function
The hard matching coefficients for the annihilation of dark-matter particles in an integer isospin-j multiplet were previously computed in [7]. For the operators O 1,2 defined in (30), (31), they read where c 2 (j) = j (j + 1) is the SU(2) Casimir of the isospin representation j, andĝ 2 (µ) denotes the SU(2) gauge coupling in the MS scheme at the scale µ. 9 They satisfy the RGE equation The one-loop anomalous dimension matrix takes the form n G = 3 is the number of SM fermion generations. For NLL resummation we also include the two-loop cusp anomalous dimension. A detailed discussion of the evolution of the Wilson coefficients after diagonalization of the anomalous dimension can be found in [7], and will not be repeated here. Details on the calculation of the hard functions and anomalous dimensions are provided in Appendix A.
It is convenient to define the vector 9 When the argument µ is omitted in the following, it is implied. Similarly forα 2 =ĝ 2 2 /(4π).
of hard functions which will be used below to demonstrate the scale invariance of the annihilation rate. The RGE for H reads with as follows from (69).

Photon jet function
The anti-collinear photon jet function is the same as for the narrow resolution case and its definition is given in (54). Since the photon jet function and the soft function have the same invariant mass squared of order m 2 W , they are defined in SCET II and require an additional rapidity regulator. We chose to use the rapidity regulator introduced in [27,28]. Details on the implementation of this regulator can be found in Appendix B.1. For completeness we report the result for Z γ ≡ Z 33 γ , already given in [7]: where ν is the scale associated with the rapidity regulator andŝ W (µ) is the sine of the weak mixing angle in the MS scheme. ∆α is the difference between the fine structure constant α = 1/137.036 and α OS ( Since Z γ depends both on µ and ν, we need to resum the photon jet function in virtuality and rapidity. We will first discuss the resummation in µ and then in ν. The RG equation is with anomalous dimension The anomalous dimensions can be expanded perturbatively in the form 10 The cusp anomalous dimension coefficients up to the two-loop order are given by with c 2 (ad) = 2 and n G = 3. The one-loop coefficient γ Zγ can be obtained from its definition. Calculating the derivative in µ of (74) using the beta-function ofŝ 2 W , which can be inferred from (323), yields In the computation of (76), we used the fact that the cusp anomalous dimension appears in the same way at all orders [29], so only a one-loop calculation is necessary to determine the prefactor of the cusp piece. Eq. (75) can easily be solved, which results in the following expression for the virtuality evolution factor where µ i and µ f denote the initial and final virtuality scales before and after evolution, respectively. Note that (80) is a general solution to the RGE (75), valid to all orders. The integral in the exponent in (80) has to be computed numerically due to the appearance of other Standard Model couplings in the β-function beyond one-loop. This is also true for the virtuality evolution factors of the other functions in the factorization theorem. More care has to be taken when performing the resummation in rapidity. The rapidity renormalization group (RRG) equation is given by with the fixed-order one-loop anomalous dimension 10 In general, starting from the two-loop order, γ i , second-order terms involving several SM couplings can appear. However, this is not the case for the cusp anomalous dimension, which is the only two-loop anomalous dimension needed at NLL'.
One could now use (82) to solve the RRG. This procedure imposes that one first evolves in rapidity and only afterwards in virtuality, because in higher orders γ ν Zγ contains terms of the form α n 2 ln m (µ/m W ) with m ≤ n. If the virtuality evolution is done first, these logarithms become large and require themselves resummation. To avoid this issue, we note that the independence of any observable of the scales µ and ν gives the condition From (75), (81) and (83) we deduce the constraint (A similar constraint also applies for the soft function, discussed in Section 3.4 below.) We can now solve (84) to obtain the integrated form of the rapidity anomalous dimension where the constant is determined such that one obtains the fixed-order non-cusp piece of the rapidity anomalous dimension, which is zero at the one-loop order (82). The logarithms ln(µ/m W ) are summed by (85) to all orders in perturbation theory. Using the integrated form (85) of the rapidity anomalous dimension, we solve the RRG (81) to obtain the resummed rapidity evolution factor where ν i and ν f denote the initial and final scales of the rapidity evolution, respectively. Expanding the argument in the exponent of V (µ, ν i , ν f ) inα 2 to order O(α 2 ), one would recover the rapidity evolution factor that can be computed from the fixed-order expression for γ ν Zγ in (82). Note that in order to confirm the µ-independence of the cross section, which will be discussed below in Section 3.5, it suffices to use this fixed-order expression. For more details on the rapidity evolution factor we refer to [28].
Depending on which resummation path is chosen, the anomalous dimensions in both evolution factors (80) and (86) are required at different order. If we first evolve in rapidity and only afterwards in virtuality, the µ-dependent logarithm in V (µ, ν i , ν f ) is never large and we only need γ cusp at the one-loop order to achieve NLL' accuracy. At the same time, the ν-dependent logarithm in U (µ i , µ f , ν) will be large and thus the virtuality evolution factor requires γ cusp at two loops. If we first resum in virtuality and then in rapidity, the situation in reversed and we need γ cusp at the two-loop order for V (µ, ν i , ν f ) and at one-loop for U (µ i , µ f , ν).
Using the resummed expression for V (µ, ν i , ν f ) and keeping in mind which order of the anomalous dimensions needs to be included ensures path independence for the µ − ν resummation, which implies the relation For the resummation of the photon jet function, we chose to resum first in rapidity and then in virtuality. As discussed above, for NLL' accuracy, this requires γ cusp in the one-loop approximation for V (µ, ν i , ν f ) and we use (86) in the form The virtuality evolution factor U (µ i , µ f , ν) is computed with the two-loop cusp and one-loop non-cusp anomalous dimension from (80). The resummed photon jet function reads Hence, the rapidity scale appearing in the virtuality RGE (80) is to be understood as the endpoint ν f ∼ m W of the rapidity evolution.

Jet function for intermediate resolution
The jet function (55) in the intermediate energy resolution case describes the unobserved hard-collinear final state with virtuality m χ m W m 2 W . It is therefore justified to neglect the masses of the electroweak gauge bosons, the fermions, and the Higgs boson, and to calculate the jet function in the unbroken regime of the SU(2) L ×U(1) Y gauge symmetry. This implies that no additional rapidity regulator is needed (contrary to the case of the narrow resolution jet function, which is further discussed in Appendices B.1 and B.2). We separately give the results of the Wilson line contribution and of the self-energy contribution, in order to better identify the origin of the different terms, and hence write The one-loop results for the unrenormalized jet function terms read where T F = T s = 1/2 and n s = 1. The jet function follows after taking the imaginary part and expanding in terms of star distributions [30]. We obtain, using (57) and the numerical values of the group factors, The definition of the star distributions is provided in (225) of Appendix B.2. The further treatment is very similar to the gluon jet function in QCD [31]. It is convenient to work with the Laplace-transformed jet function j int since it renormalizes multiplicatively. The Laplace transform of J int (p 2 , µ) is defined by where l = 1/(e γ E τ 2 ) and the explicit result after renormalization reads The corresponding RG equation is the ordinary differential equation with Laplace-space anomalous dimension γ J is needed at the one-loop order for NLL' resummation, The RGE (96) is solved by where µ j ∼ √ m χ m W is the natural scale of the hard-collinear jet function and the integrals S(µ j , µ) and A γ J (µ j , µ) are defined as The variable η is defined by As mentioned before, at NLL' the integrals S(µ j , µ) and A γ J (µ j , µ) can only be solved numerically due to the appearence of several SM couplings in the β-function for α 2 beyond one loop. Note that in the second line of (99), the logarithm in the argument of the Laplace-transformed jet function has been traded for a derivative with respect to η. The complete τ -dependence of j int is then contained in the factor (τ 2 /µ 2 j ) η , and the inverse Laplace transform becomes simple. By exploiting the relation one obtains the resummed jet function in momentum space in the form The more complicated jet function for the narrow energy resolution, J nrw , was used in [7]. In Appendix B.2, we provide details on its computation and give its full expression.

Soft function
The soft function is defined by the vacuum amplitude (59) of the soft operator (53) with index contraction as specified in (60). Let us recall that the soft operator is the product of soft Wilson lines arising from the decoupling of soft SU(2) gauge bosons from the four particles in the 2 → 2 annihilation amplitude. The SU(2) indices are then contracted in a way that depends on the operator O i and the external DM two-particle state I, resulting in the function W ij IJ (ω). The soft function is sensitive to physics at virtualities of order m 2 W , and therefore must be computed in the effective theory with broken SU(2) L ×U(1) Y gauge symmetry and massive SM particles (unless the mass is much smaller than m W ). For NLL' accuracy, the one-loop soft function and its NLL RG evolution is needed. Here we summarize these results. Technical details on the computation of virtual and real one-loop diagrams are given in Appendix C, including the regularization, which involves rapidity regularization, together with some observations on partial virtual-real singularity cancellations.
The virtual one-loop contributions to the soft function are the same as for the narrow resolution case and they were already computed in [7]. In the intermediate resolution range, the real emission of soft EW gauge bosons is kinematically allowed. The new contributions as well as the virtual diagram results are given explicitly in Appendix C.
To guide the discussion we present here the result for the W 22 (+−)(+−) component of the soft function, which has the most complicated structure and allows us to explain the resummation procedure: The complete set of soft function components for all operator and two-particle-state combinations is collected in Appendix C.4. As for the unobserved-jet function, renormalization becomes multiplicative in Laplace space. The forward and inverse Laplace transforms are defined as As can be seen from (105), the Laplace transforms required for the soft function are where the functions si, ci are defined as It is convenient to introduce the following vector notation for the Laplace transformed soft functions. The RRG equations for the soft functions take the form where the fixed-order one-loop rapidity anomalous dimension is given by Note that the non-cusp piece of Γ ν W is zero at one loop. The discussion of the rapidity evolution factor from Section 3.2 equally applies to the soft function. We hence use (84), (85) and (111) to compute the rapidity-resummed soft function where Γ ν W (µ) is the integrated rapidity anomalous dimension for the soft function. As was discussed in the case of the photon jet function, the order of the anomalous dimensions included in (113) depends on the resummation path in the µ − ν plane. The RRG is diagonal in both the operator index, encapsulated in the vector notation, and the twoparticle state index pair IJ. Notice that only the soft functions and the photon jet function depend on the rapidity scale. We make the choice to evolve the photon jet function from the jet rapidity scale ν h ∼ 2m χ down to ν s ∼ m W . This means that we can set ν = ν s for the soft function, which makes the rapidity evolution factor (113) equal unity.
The virtuality RG equation for the Laplace-transformed soft function is also diagonal in IJ, but its non-cusp piece is non-diagonal in operator space, with anomalous dimension As in the case of the photon jet function, (112) and (115) can be obtained from their definitions by taking the derivatives in µ and ν, respectively, of w IJ . At the one-loop order, which is enough for NLL' resummation, the anomalous dimension γ W evaluates to γ (0) The solution to (114) takes the form The evolution matrix U W is diagonal, and the diagonalization matrix R and its inverse R −1 are given by The integrals A γ W and η have been introduced in (101) and as already explained there, they can only be solved numerically at NLL'. As a last step we need to go back to momentum space and compute the inverse Laplace transform of (117). The entire dependence on κ is contained in w IJ (s, µ s , ∂ η ) κ ν η . We therefore defineˆ W IJ (ω, µ s , ν) to be the inverse Laplace transform of w IJ (s, µ s , ∂ η ) κ ν η : The inverse transformation requires the computation of For the above representative index and operator combination IJ = (+−)(+−) and ij = 22, the inverse Laplace transform giveŝ The results forŴ ij IJ in all possible index and operator combinations IJ and ij are collected in Appendix C.5. Finally, using (117) and (120), we find that the virtuality resummed soft function in momentum space takes the form We emphasize that we did not include the rapidity evolution factor (113) in (124), since we evolve Z γ in ν from ν h to ν s which makes the soft function rapidity evolution factor unity.

RG and RRG invariance of the cross section
The factorization formula for the intermediate resolution case given in (40) and (41) puts constraints on the anomalous dimensions, since the physical photon energy spectrum has to be independent of the virtuality and rapidity factorization scales µ and ν. This independence on the scales manifests itself in the two consistency equations Note that the Sommerfeld factor (17) is computed at leading order, which makes it scale independent so it does not have to be taken into account when computing (125) and (126). In previous subsections, we already made use of the fact that a Laplace transformation turns convolution into multiplication. It is thus easiest to derive the implications of (125) and (126) in Laplace space, by taking the Laplace transform of (40), (41) with respect to the variable e γ ≡ 2(m χ − E γ ). Calling the Laplace variable t, the Laplace transform of the convolution of the jet with the soft function is (for brevity, we omit µ and ν in the arguments, as well as operator and two-particle state indices) When going from the second to the third line in (127), we made use of the substitution p 2 = 2m χ (e γ − ω). Also, since p 2 is strictly positive, we can set the lower p 2 -integration boundary to zero. Using the definitions (94), (106) of the Laplace-transformed jet function and soft function, respectively, we arrive at the fourth line of (127). We can therefore write (125) as Taking the derivative and making use of the definitions (72), (75), (96) and (114) of the anomalous dimensions results in The terms in (129)   Left: common reference scale is µ s . Right: common reference scale is µ j . In both cases, Z γ is evolved in ν from ν h to ν s .
The same steps can be applied for the evaluation of (126), except that in (128) we differentiate with respect to ln ν. Since only the photon jet function and the soft function depend on the rapidity scale ν, using the definitions (81) and (111) results in the rapidity consistency equation This can be shown to be satisfied by the values for the rapidity anomalous dimensions given in (82) and (112).
Since (129) and (131) are fulfilled, we confirm that at the one-loop order the factorized cross section is independent of the scales µ and ν. It should be noted that the cancellation of the off-diagonal non-cusp terms of Γ H and Γ µ W in (129) is non-trivial. In total, this provides a strong check of the consistency of the calculation. The corresponding consistency check for the factorization of the narrow resolution case is presented in Appendix D.

Resummation schemes
Having collected the RG equations for all the factors in the factorization formula, we show in Figure 2 two different possibilities for the resummation of the functions appearing in the factorization theorem. For the first resummation scheme, shown in the left Figure 2, we choose µ s as the common reference scale and evolve the Wilson coefficients C i and the unobserved-jet function J down to the soft scale µ s , while the soft function W and the photon jet function Z γ do not contain large logarithms when evaluated with µ = µ s , and hence do not require resummation in µ. Resummation in the rapidity scale is however necessary. We choose to evolve the photon jet function from ν h to ν s . Equivalently one could also evolve the soft function from ν s to ν h . This resummation scheme is close to the implementation of the narrow resolution case [7], where there is no hard-collinear scale µ j , and the hard functions are evolved all the way from the hard to the soft scale.
A more conventional implementation of resummation in the presence of an intermediate hard-collinear scale is the second resummation scheme illustrated in the right Figure 2. Here we choose µ j as the common reference scale, and evolve C i down, and Z γ and W up to µ j . Z γ is evolved in rapidity from ν h to ν s as before. Note that in this second case, as discussed in Section 3.2, the specific form of the rapidity evolution factor V depends on whether we first evolve in ν and then in µ or vice versa. Since we saw that V takes a simpler form if we resum first in ν and then in µ, we choose this ordering, as is also shown in Figure 2 (right).
Both schemes give the same results up to effects beyond the accuracy of the truncation of the RG equations.

Results
In this section we present the results for the DM annihilation process χ 0 χ 0 → γ + X, assuming an intermediate energy resolution E γ res of the instrument of order of the weak scale m W . First we show σv (E γ res ), as defined in (1), as a function of the DM mass m χ and then perform a numerical comparison of the present calculation with the narrowresolution result of [7]. An analytic comparison of the two energy resolution cases is made in the next section, where we discuss the logarithms in the annihilation rates up to the two-loop order.

Energy spectrum
The upper panel of Figure 3 shows the cumulative endpoint annihilation rate σv (E γ res ), plotted as a function of the DM mass m χ . The mass range includes the first two Sommerfeld resonances. The different lines refer to: the Sommerfeld-only calculation (blackdotted), also called "tree", since Γ IJ is evaluated in the tree approximation without any resummation, and multiplied with the Sommerfeld factor S IJ according to (40); the LL (magenta-dotted-dashed), the NLL (blue-dashed) and finally the NLL' (red-solid) resummed expression for Γ IJ , the latter of which represents the calculation with the highest accuracy. The photon energy resolution is set to E γ res = m W in this figure. The lower panel of the Figure shows the same LL, NLL and NLL' resummed annihilation rates, but normalized to the Sommerfeld-only result for better visibility of the resummation effect. We see that the resummation leads to a substantial reduction of the cross section, as is generally expected for Sudakov resummation. The size of the effect is consistent with the finding of previous computations [3][4][5]7] of related observables or in different resolution regimes. In particular, in the interesting mass range around 3 TeV where wino DM accounts for the observed relic density, the rate is suppressed by about 30 − 40%.
The resummed predictions are shown with theoretical uncertainty bands computed from a parameter scan with simultaneous variations of all scales. Specifically, the scales µ h , ν h were varied in the interval 2m χ [1/2, 2], µ j was varied in the interval 2m χ m W [1/2, 2] and µ s , ν s were varied in the interval m W [1/2, 2]. The errors were then determined very conservatively by taking the maximum and minimum values in this five-dimensional parameter space. This scan was repeated for each mass point. For each parameter scan, we specified 21 values distributed logarithmically in the intervals given above, with ten values above and ten below the central values of the intervals.
We find that the residual theoretical uncertainty at the NLL' order becomes negligible and is given by the width of the red-solid curve in Figure 3. It is also apparent that the different levels of resummation successively reduce the theoretical uncertainty considerably, from 15% at LL, to 9% at NLL and 1% at NLL' at m χ = 2 TeV. Numerically, for the two mass values m χ = 2 TeV (10 TeV) the ratio to the Sommerfeld-only rate is 0.641 +0.103 It is instructive to separate the integrated photon energy spectrum σv (E γ res ) into the contributions due to the different Sommerfeld factors in (40). Thus, we write for the smaller DM mass value m χ = 500 GeV. We observe that at m χ = 2 TeV (and similarly for larger masses), the Sommerfeld factors are large, as expected, and the annihilation rate is dominated by the (+−)(+−) hard annihilation channel, which starts at tree level in the fixed-order expansion. The Sommerfeld factors are O(1) and even smaller than 1 for m χ = 500 GeV for the off-diagonal annihilation contributions (00)(+−) and (+−)(+−), for which the Sommerfeld enhancement does not compensate the loop suppression at small masses. The results shown in this section were computed with the more conventional second of the two resummation schemes discussed in Section 3.6. We implemented both schemes and found full numerical agreement at NLL' at the 0.1% level, as also follows from the analytic comparison, see (160) below.

Matching energy resolutions
In the introduction we identified three different regimes for the energy resolution E γ res , the narrow, the intermediate and the wide region. These cover the entire range of E γ res for DM indirect detection experiments. In [7] we provided NLL' predictions for the photon-energy spectrum near the endpoint assuming a narrow energy resolution of E γ res ∼ m 2 W /m χ , close to the line signal, while in this work we focus on E γ res ∼ m W , which is more realistic for present and future indirect DM searches in the TeV energy region. The two calculations differ in the structure of the unobserved jet function and the soft function, and exhibit different large logarithms. The question arises whether the two computations can be matched to provide an accurate result for the entire range from E γ res ∼ 0 to E γ res ≈ 4m W , which we tentatively define as the upper limit of validity of the intermediate resolution case.
In Figure 4, we show the annihilation cross sections for the narrow (blue-dotted) and the intermediate resolution We observe a wide interval in E γ res , covering the range of resolution in between the validity regions of the two calculations, for which the annihilation rates in both calculations agree with high precision. At low resolution there is a steep rise of the narrow resolution rate, which occurs at E γ res ≈ 4m 2 Z /m χ . Above this value the resolution is not enough to separate the γZ contribution, leading to a sharp increase of the semi-inclusive rate. Since the unobserved-jet function for the intermediate resolution cross section is computed under the assumption that the particles are massless, this feature is absent in this curve (dashed/red), which is hence clearly not valid for very small resolution. In the narrow resolution regime the invariant mass of the unobserved-jet function also passes through the W + W − , ZH and tt thresholds. However, these thresholds are not visible on the scale of the plot. The narrow resolution computation agrees very well with the inter- mediate one well into the regime of validity of the latter, and vice versa. As one moves to even higher E γ res , the intermediate resolution line starts to depart from the narrow resolution one. Here the narrow resolution computation clearly ceases to be accurate, because it fails to capture the effect of soft electroweak gauge boson radiation, which is now kinematically allowed. Nevertheless, even at the highest E γ res = 1 TeV shown in the plot, the difference stays below the 20% level, as can be seen from the ratio plots at the bottom of the two panels in Figure 4. We note that as the DM mass becomes larger, the separation between the two validity regimes (the shaded bands in the Figure) increases, but the matching continues to work well even for the 10 TeV DM mass example.
These observations show that the present work and [7] combined result in highly accurate theoretical predictions for the photon energy spectrum in dark matter annihilation, here for the wino model, in the entire energy resolution range from E γ res ∼ 0 to E γ res ≈ 4m W . It would be interesting to perform a similar matching between the results of the present paper and the results of [6], which would extend the knowledge of the resummed energy spectrum to even wider resolution. As discussed in the introduction, with the anticipated energy resolution of the CTA experiment, we expect this to be necessary for DM searches only in the 10 TeV mass region and beyond.

Fixed-order expansions
In this section we perform analytic expansions of the annihilation rate matrix Γ IJ up to the two-loop order. This provides some insight into the structure of large logarithms in the photon energy spectrum at large photon energy, depending on the energy resolution, and explains why the two computations agree remarkably well over a large interval of E γ res , as observed in the previous section. Readers interested only in the numerical result for the spectrum may skip this section.

Double-logarithmic approximation
Before moving to fixed-order expansions it is instructive to compare the NLL' result to the double-logarithmic approximation. This approximation is obtained by a) evaluating all functions in the tree-approximation, b) keeping only theα 2 × log 2 terms in the exponents of the RG evolution factors. For the two resolution regimes discussed in this paper, the double-logarithmic approximations read σv int (E γ res ) = The dependence of the coefficient of large logarithms on the energy resolution is already apparent from these equations. Since the 'nrw' formula describes an observable that , which partially compensates the Sudakov suppression associated with the hard-function resummation.
The double-logarithmic approximation is visualized in Figure 5. It is seen that within their respective validity ranges (shaded areas in the plots) the double-log approximations of the intermediate and narrow resolution results are close to the full NLL' resummed results, shown for comparison (dimmer dashed/red and dotted/blue curves). In the narrow resolution case the step function in (136) correctly describes the sharp rise of the annihilation cross section due to the opening of the γZ channel. However, Figure 5 also demonstrates that the precise shape of the cumulative annihilation rate in E γ res and, in particular, the smooth matching of the two resolution regimes observed in the previous section cannot be explained in the double-logarithmic approximations (136), (137). We therefore analyze the subleading logarithms in the oneand two-loop order in the following subsection.

Expansion of the resummed annihilation rate
We re-expand the resummed annihilation rates (41) where, by construction, the coefficients c (n,m) IJ (E γ res , µ) are O(1) numbers, and the large logarithms ln(2m χ /m W ) are made explicit. Note that the coefficients c (n,m) IJ (E γ res , µ) are different for the two resolution regimes. 11 The resummed rate depends on many scales, µ from the renormalization of the coupling, and the scales from the initial and final values of the RG evolution. To make the large logarithms explicit, we normalize scales by their natural values. For example, ln(µ 2 j /m 2 W ) is written as ln(µ 2 j /(2m χ m W ))+ 1 2 ln(4m 2 χ /m 2 W ), such that the first logarithm is O(1) and part of one of the c Before discussing the behaviour of the c (n,m) IJ coefficients as functions of E γ res , let us clarify which logarithms in (138) are captured by NLL' resummation. After RG evolution, the resummed annihilation cross sections are obtained in the form with functions f i (α 2 L) of the O(1) quantityα 2 L ≡α 2 ln(4m 2 χ /m 2 W ). The LL approximation amounts to keeping f 0 , NLL adds f 1 , while NLL' adds C 1 . Other terms not 11 In the following we drop the arguments E γ res , µ for brevity. 12 We use the following terminology: "n-loop" refers to the O(α n 2 ) correction only, while NLO refers to the sum of tree and one-loop, etc. written are beyond the NLL' accuracy. Expanding inα 2 , we observe that NLL' resummation determines the three highest powers of logarithms in any order of perturbation theory, specifically c in (138) for all n. In particular, for n = 1 (one-loop) the NLL' resummation determines all the possible coefficients that exist at this order, including the non-logarithmic term m = 0, while at two loops (n = 2) all logarithms except the single logarithm are obtained. Since the dependence on the matching scales such as µ j introduced by resummation must cancel at every fixed order, those fixed-order coefficients, which are obtained exactly from expanding the resummation formula, must be independent of these scales. On the other hand, at two loops, the single logarithmic and constant terms still depend on O(1) quantities such as ln(µ 2 j /(2m χ m W )) as can be seen from the explicit expressions in Appendix E.
In the following we discuss the logarithmic structure for the channel IJ = (+−)(+−), which is the most interesting one, since the other channels do not have a tree-level coefficient. 13 We then evaluate the coefficients outside of their validity range, for example we take a coefficient from the double Taylor expansion of the narrow resolution formula and extrapolate it to E γ res ∼ m W or to the transition energy resolution scale E γ res ∼ (m W /m χ ) 1/2 m W in order to study the numerical matching of the two resolution cases. The extrapolation induces a reshuffling of the logarithms in (138) because O(1) coefficients in one regime may develop large logarithms in the other.

Tree level
The tree-level coefficients in (138) are c c int(0,0) The narrow resolution formula distinguishes the contribution from the γZ line from the γγ while the intermediate resolution formula does not. When evaluated at E γ res > m 2 Z /(4m χ ) both formulas yield the same result.

One loop
The one-loop term in (138) reads explicitly For the presentation of the coefficients, some abbreviations will be helpful: The variable L is the large logarithm in the expansion (138). Note that l R is an O(1) quantity for intermediate resolution, but counts as a large logarithm in the narrow resolution case. The fixed-order expansion is performed in the running couplingsα 2 (µ), s 2 W (µ) at the scale µ of order m W . We define the O(1) quantity l µ 2 ≡ ln(µ 2 /m 2 W ). These explicit µ-dependent logarithms cancel the implicit scale dependence of the couplings up to residual dependence of higher order than the NLL' accuracy of the approximation. In addition to the variables introduced above, we define and the resolution-dependent function j(E γ res ) by means of the equation where the one-loop contributions to Z γ (µ, ν) and J 33 nrw (p 2 , µ, ν) are given in (74) and (220), respectively. The function j(E γ res ) captures the complicated dependence of J 33 nrw (p 2 ) on the masses of the SM particles and E γ res (see Appendix B.2), and is constructed such that it is independent of µ and ν.
With these abbreviations at hand, we find [σv] int 1−loop where c nrw(1,0) The l µ 2 dependence in the c (1,0) (+−)(+−) coefficients above is compensated by the running couplingsα 2 (µ) andŝ 2 W (µ) in the corresponding LO terms. The coefficients depend on E γ res through the functions l R , λ R and j defined above. Therefore, in order to investigate the transition from the narrow to the intermediate resolution formulas we need to understand the asymptotic behaviour of these functions. For instance, up to corrections of order m 2 W /(4m χ E γ res ). This can be obtained from expanding the explicit expressions for the one-loop Wilson line and self-energy contributions given in Appendix B.2, or, more simply, by performing the expansion by regions [16] before taking the integrals.
When extrapolating (145) into the intermediate resolution regime, we can write j(E γ res ), using ln(4m χ E γ res /m 2 W ) = 1 2 L + l R , as Then, for E γ res m 2 W /m χ , We note that due to the asymptotic behaviour of j(E γ res ), the large logarithms precisely match. The difference is a non-logarithmic term, which turns out to be quite small, and amounts to O(1%) of the tree-level cross sections independent of the DM mass. This is visualized in Figure 7 where the one-loop coefficient (excluding the factorα 2 /π) is plotted for the two resolutions (narrow in dashed/blue, intermediate in solid/red). The absolute value of these dimensionless coefficients is, for both cases, large but the coefficients differ by no more than 3% in the hatched cross-over region. Similar results are found when I or J = (00) as can be verified from the coefficients c (n,m) IJ listed in Appendix E. One may wonder why in (151) the narrow and intermediate resolution coefficients do not agree exactly, since by construction the NLL' approximation reproduces the full one-loop calculation. However, this is true only up to power corrections in m W /m χ . The difference in (151) arises from the λ R term in the intermediate resolution coefficient (146). In the narrow resolution limit λ R is a power-suppressed effect of order m W /m χ .

Two loops
The two-loop term in (138) reads c int(2,2) for the intermediate resolution case. As before, the l µ 2 dependence of the c  In order to understand the behaviour of the curves in the Figure analytically, we use the asymptotic behaviour of the E γ res -dependent functions within the coefficients. The leading-logarithmic dependence of these for different E γ res scaling is shown in Table 1. Besides the two energy resolution regimes associated with the 'nrw' and 'int' factorization theorems, the transition scale constructed from the geometric mean of the narrow and intermediate resolution scales is also considered.
In Table 2 we show the leading logarithm that results from reevaluating the 'nrw' and 'int' two-loop coefficients at the three scales in E γ res relevant to Figure 8. We verify the behaviour encountered in both panels of the Figure. Namely, in Figure 8 the two-loop Table 1: Leading-logarithmic dependence of the E γ res -dependent functions appearing in the fixed-order expansions when evaluated at the three E γ res -scales relevant to Figure 8. Vanishing entries are to be understood as power-suppressed. coefficient associated with the 'nrw' formula is larger than the corresponding one from the 'int' formula for the narrow resolution regime. This property is supported by the positive difference of the leading logarithmic term (L 4 /32) between the two formulas as evaluated at the 'nrw' regime (last row and second column of Table 2). Conversely, when E γ res ∼ m W the opposite happens, and the 'int' coefficient is larger than the 'nrw' one, consistent with the last entry (last column from left to right and last row from the top to the bottom) of Table 2. The vanishing of the O(L 4 ) term for E γ res ∼ m W m W /m χ explains the almost perfect matching of the 'nrw' and 'int' coefficients in the transition region as observed in Figure 8.
In summary, that the matching works so well over a wide range of energy resolution is a consequence of the smallness of the difference in the leading logarithms. For example, extrapolating the narrow resolution coefficient to intermediate resolution, we find at NNLO At one loop, as discussed before, the difference lacks large logarithms since λ R is an O(1) function of E γ res provided E γ res ∼ m W . At the two-loop level we see a partial cancellation of the L 4 coefficients (as 1/32 1). In (159) we therefore include the L 3 term, which constitutes the largest difference term at two loops when L is not extremely large. Table 2: Leading-logarithmic terms of the two-loop coefficients in (138) for the 'nrw' and 'int' factorization formulas, and the difference of the two, at the scales relevant to Figure 8.

Resummation schemes compared
So far the discussion on the fixed-order expansions of the intermediate resolution formula has been done using the first resummation scheme of Section 3.6. This is the most natural choice when comparing with the factorization formula in the narrow resolution case. We performed the same fixed-order analysis for the second resummation scheme and found exact agreement in all the coefficients at two loops except c [σv] tree where ϕ R is defined in (310). Numerically this difference is not be larger than O(0.1%) of the tree-level cross section. Note that since the single log coefficient is not obtained unambiguously by NLL' resummation, also c (2,1) (+−)(+−) could have depended on the resummation scheme, but this turns out not to be the case.

Conclusion
The search for high-energy photons plays an important role in detecting dark matter through its annihilation in the center of the Milky Way, or in dwarf galaxies. Connecting a possible signal to a DM model, or to place limits on the parameters of the model, including the DM mass itself, requires an accurate theoretical calculation of the annihilation rate. When the DM particle carries electroweak charges and its mass is much larger than the mass of the electroweak gauge bosons, standard perturbation theory in the small couplings of the SM breaks down. Large enhancements of loop diagrams due to non-relativistic scattering and due to soft and collinear gauge bosons must be summed to all orders in the coupling expansion. In this paper we considered the photon energy spectrum of the semi-inclusive photon final state γ + X, integrated from the endpoint E γ = m χ over an interval of size E γ res , which corresponds to the observable measured by γ-telescopes, when the flat integration in (1) is replaced by the instrument-specific resolution function of characteristic width E γ res . The main theoretical result is the factorization formula (40) for the annihilation rate for energy resolution E γ res ∼ m W (41), and E γ res ∼ m 2 W /m χ (42), respectively, and the calculation of the all-order resummed rate to NLL' accuracy in the electroweak Sudakov logarithms. The main results relevant to observations are summarized in Figures 3 and 4. The corresponding result for narrow resolution has already been shown in [7], but the derivation and the matching to the intermediate energy resolution E γ res ∼ m W has been presented here for the first time. While the theoretical formalism is more general, and so are some of the calculations, the complete NLL' calculation has been performed in the so-called pure wino model, where the SM is extended by a fermionic SU(2) triplet, of which the electrically neutral member is the DM particle. We highlight the following two observations: • Electroweak Sudakov effects are large and reduce the annihilation rate to highenergy photons by about a factor of two in the multi-TeV region. As soon as the full one-loop effects are included, that is, the accuracy of the calculation elevated from NLL to NLL', the theoretical uncertainty, as measured by renormalization and factorization scale variation, becomes negligible (about or below 1%), see Figure 3.
• The two separate calculations for narrow and intermediate energy resolution match very accurately, resulting in precise theoretical results from the line-like final state at E γ res ∼ 0 to E γ res ∼ 4m W (perhaps, beyond), see Figure 4. While the calculations apply to any DM mass with m χ m W , given the energy resolution of the H.E.S.S. and CTA experiments, they are most relevant for m χ in the range between 1 and 10 TeV. This is also the range where the wino model is most compelling.
In [6,8] a complementary approach has been pursued, which applies to what we called "wide" energy resolution E γ res m W . The available results are of NLL accuracy for the same wino model, and, given the observations above, it would be of interest to a) extend them to NLL' and b) match them to the intermediate resolution case discussed here.
The results shown here demonstrate the success of EFT techniques, non-relativistic and soft-collinear, to deal with the breakdown of electroweak perturbation theory in the high-energy regime. This opens the perspective to extend the calculations to models other than the wino model. Given the small uncertainty of ≤ 1% from scale variation of the resummed perturbative expansion, it is probable that the largest theoretical uncertainty now arises from modifications of the Sommerfeld effect due to sub-leading effects in the non-relativistic effective theory, and, for smaller m χ , from power-suppressed effects of order m W /m χ , which are systematically neglected in the present treatment.

Acknowledgement
This work has been supported in part by the Collaborative Research Center 'Neutrinos and Dark Matter in Astro-and Particle Physics' (SFB 1258) and the Excellence Cluster 'Universe' of the Deutsche Forschungsgemeinschaft. AB acknowledges the support by the ERC Starting Grant REINVENT-714788.

A Hard matching coefficients
In this Appendix we provide more details on the calculation of the hard-matching coefficients given in (67) and (68).

A.1 Amplitude in the full theory
The matching condition between the full theory (SM plus an isopsin-j dark matter multiplet) and the effective theory requires that the on-shell amplitudes for 2 → 2 annihilation of two dark matter fields to two SU(2) gauge bosons computed in the two theories must be equal: Here the left-hand side refers to the UV-renormalized amplitude in the full theory. The symbol {p, s} indicates the dependence of the amplitudes on the momenta and the spin/polarization orientations of the four external particles. The operators O i are Swave operators. To extract their coefficient we can set the relative momenta of the annihilating particles to zero. We choose p 1 = p 2 = m χ (1, 0) for the initial state, and p 3 = m χ n − , p 4 = m χ n + for the final state. We define projectors applied to the full theory amplitude such that where M i, full are the full-theory projected amplitudes corresponding to the gauge and spin structures of the two operators O 1,2 defined in (30) and (31). The expressions in (162) directly correspond to the bare matching coefficients since the loop diagrams in the effective theory are all scaleless and vanish in dimensional regularization. The two projectors have the explicit expressions where = (4 − d)/2, d is the space-time dimension and c 2 (j) = j(j + 1) for an isospinj representation. The projectors differ only in the SU(2) part, since both operators have the same Dirac and Lorentz index structure which projects only on the spin-singlet contribution of the amplitude. The projectors can be considered as operators in spin space and the same is true for the amplitude. We compute the matching coefficients at the one-loop order (see Figure 9 for a sample of diagrams). We use dimensional regularization for both ultraviolet and infrared singularities. The calculation of the bare full theory amplitudes has been carried out by using a set of computer-algebra tools. FeynRules [32], FeynArts [33] and FormCalc [34] were used in combination for the model implementation and the amplitude generation. The algebraic manipulations and simplifications have been carried out with a private code written in FORM [35]. The reduction to master integrals at threshold was performed with Reduze [36]. We calculate the Feynman diagrams in the unbroken SU(2) gauge theory and find for the bare projected full-theory amplitudes where We find the same expressions, both for the case of Dirac and Majorana fermions. This is no surprise since a possible difference could arise only from s-channel diagrams with a fermion-fermion-gauge boson vertex. At threshold these diagrams do not contribute to the amplitude. In the following, we find it convenient to suppress the µ dependence of the renormalized SU(2) coupling in intermediate results. We remove the UV divergences by coupling, field and DM mass renormalization. The coupling constant is renormalized in the MS scheme while the mass and field renormalization is done in the on-shell scheme so that no further residue factor is required to obtain the on-shell amplitude. The SU(2) coupling, DM mass and field renormalization, and the SU(2) gauge boson field renormalization constants are given, respectively, by where c(j) = c 2 (j)(2j + 1)/3 and n G = 3 is the number of fermion generations. In (166) the term 2c(j)r/3 corresponds to the heavy DM fermion contribution, the term −43/12 to the gauge boson and Higgs contributions, while the 2n G /3 piece arises from the SM fermion loops. The parameter r assumes the values r = {1, 1/2} for Dirac and Majorana fermions, respectively. In the effective theory the heavy fermion is integrated out and does not contribute to the running of the gauge coupling anymore. Hence, similarly to switching between schemes with different massless quark flavours in QCD, we decouple the DM contribution from the running of the gauge couplingĝ 2 through the substitution in (164). After this replacement the r dependence drops out, and the final result will be independent of the Dirac or Majorana nature of the fermion. The UV-renormalized projected full-theory amplitudes, which equal the bare Wilson coefficients, read The remaining IR divergences must be cancelled by matching. This will be done in the next subsection by operator renormalization in the effective theory.

A.2 Operator renormalization in the effective theory
In the effective theory the loop diagrams are all scaleless and therefore vanish in dimensional regularization. The EFT matrix element is therefore given by the tree matrix element and the tree diagrams with counterterm insertions. To compute the UV counterterms in the effective theory, we need to regulate the IR-divergences with a different regulator than dimensional regularization. To this purpose we take slightly off-shell momenta for the incoming DM fermions and the final-state gauge bosons. By direct calculation of the effective theory diagrams we obtain for the UV poles in the MS scheme where O i tree correspond to the tree-level matrix elements of the first or second operator. Notice that the divergent parts shown do not depend on the infrared regulator and that they only depend on the hard scale 2m χ . We still need to add the external field MS renormalization factors for the effective theory fields, which read By combining everything we arrive at Coupling renormalization contributes only at higher orders inĝ 2 . The MS operator renormalization constants Z ij are a matrix in operator space such thatÔ bare i = Z ijÔ ren j (µ), i, j = 1, 2. We obtain from the above By making use of the matching condition in (161), the decoupling relation in (170) and we find that all 1/ poles cancel and we obtain the explicit results for the hard matching coefficients given in (67) and (68).

A.3 Operator Z-factors from the anomalous dimension
A second way to obtain the operator renormalization Z ij factor is to adapt the anomalous dimension known for QCD processes [37,38] to the SU(2) gauge group. We switch to the operator basis where the DM bilinear is in a definite isospin representation (the DM bilinear can be either in a singlet or in a quintuplet representation) The advantage of this basis is that the anomalous dimension at threshold is diagonal [38], where c 2 (ad) is the Casimir value of the gauge boson in the adjoint representation, and c 2 (J) the one for the DM fermion pair in the representation J = 0 (singlet) or J = 2 (quintuplet). The quantity γ ad is the gauge boson anomalous dimension and γ J H,s is the anomalous dimension of the heavy fermion pair. The anomalous dimensions have perturbative expansions in terms ofα 2 (and, possibly, other couplings in higher orders than given) The Higgs contribution −16/9 to the two-loop cusp anomalous dimension has been extracted from the -scalar contribution computed in [39].
The operator Z-factor in the MS scheme can be obtained from the anomalous dimension. Up to orderĝ 2 2 it reads where Γ = −2c 2 (j)γ cusp . In the diagonal basis defined in (180) we find Transforming back to the unprimed basis (30), (31) we find agreement with (178).

B Collinear functions and rapidity regularization
In this Appendix provide some details on the rapidity regularization, which is required for collinear and soft functions, on collinear integrals, and we supply the lengthy expressions for the narrow resolution jet function that were not given in [7].

B.1 Rapidity regularization
We employ the rapidity regulator introduced in [28], which amounts to the following replacements in the eikonal Feynman rules that originate from soft and (anti-) collinear Wilson lines collinear emission : anti-collinear emission : soft emission from (anti-) collinear direction : soft emission from the heavy line : η is the rapidity regulator and ν is a newly introduced rapidity scale, the equivalent of µ in dimensional regularization. Notice that the rapidity regulator in the above expressions is consistent between soft and collinear integrals since 2k 3 → n + k (2k 3 → n − k) in the (anti-) collinear limit. In our case all the soft and photon jet functions always require rapidity regularization, but the unobserved-jet function only in the narrow resolution case. In the following we focus on the collinear and anti-collinear scalar integrals, which appear, respectively, in the unobserved jet and photon jet function calculation. In Appendix C we present the computation of the relevant soft virtual and real integrals.
As an example we compute the off-shell collinear scalar integral (p 2 = 0) relevant to the unobserved-jet function in the narrow resolution regime. From contour integration we find that n + k < 0, such that we can replace the absolute value by −n + k.
We proceed with the loop integration by introducing the Feynman parametrization After performing the loop integration we arrive at (196) We make the substitution x 1 → x 1 x 2 and integrate first over x 2 . For convenience we drop the +iε in the intermediate expressions, and we identify p 2 → p 2 + iε as follows from the definition of the integral. After integrating over x 2 we obtain We rewrite part of the integrand as where r ≡ m 2 W /(−p 2 ), and perform the variable change x 1 = (1−y)/y. The y integration amounts to where F 1 is the F 1 -Appell hypergeometric function. This gives I c (p 2 ) to all orders in η and . Numerical checks were performed to ensure the consistency of (196) and (199). We need to expand the result first in η → 0, using the formula and then in . The +-distributions acting on a test function are defined as For the y-integral (199) we find We find this compact result after introducing the variables and with the help of relations between polylogarithms of different arguments in intermediate steps. We used the package NumExp [40] for the numerical expansion of the Appell F 1 function in (199) to check (202). In total we find The integral above is real in the region p 2 < 0, but we need to extract potential imaginary parts in the regions p 2 > 4m 2 W and 0 < p 2 < 4m 2 W . To obtain the result in the region p 2 > 4m 2 W from (204) one needs to perform the substitution In the region 0 < p 2 < 4m 2 W the result does not develop an imaginary part and it can be obtained from (204) by making the substitution where we defineβ In the intermediate resolution case the external momentum p has hard-collinear scaling p µ ∼ m χ (λ, 1, √ λ) such that p 2 ∼ λm 2 χ , while the square of the gauge boson mass scales as m 2 W ∼ λ 2 m 2 χ . Hence the expansion for p 2 m 2 W becomes relevant. Directly expanding the result in (204) yields up to power corrections, which seems to be at variance with the gauge-boson mass independence of the result for the hard-collinear jet function in the main text. However, the integral (194) is now a two-scale object. We find that there are two regions contributing to this integral, namely the hard-collinear and the soft region, k µ ∼ (λ, λ, λ). To extract the soft contribution, we need to expand the propagator in this region. The soft integral then reads Calculating the integral in a similar way as above, we obtain The rapidity regulator is only needed in the soft contribution and not in the hard-collinear one. In the hard-collinear region we can drop the gauge boson mass at leading power, and the integral evaluates to By adding up the two contributions (211) and (212) we reproduce (208). After dressing the collinear-soft scalar integral (210) with the proper tree-level factors to obtain the soft contribution to the jet function Wilson line diagram, and after taking its imaginary part, we find that the virtual (single-particle cut) piece evaluates to a scaleless integral while the real emission (two-particle cut) piece is non-vanishing. It can be shown by direct comparison that this last term equals the soft emission diagram in (269) after the convolution with the tree-level jet function has been done. This shows that in the small mass limit, m 2 W p 2 , the soft region of the jet function integral is correctly reproduced by the soft function in the factorization formula for intermediate resolution and should not be assigned to a mass-dependent collinear function. The hard-collinear region only contributes the mass-independent unobserved-jet function.
The photon jet function also requires rapidity regularization. In this case only virtual diagrams contribute. As an example, we compute the rapidity and dimensionally regulated on-shell anti-collinear scalar integral which can be obtained from (194) by setting p 2 = 0 and replacing n + k → n − k. We parametrize the integration measure by and rewrite the integrand as We perform the n + k integral first by closing the contour in the upper half plane and pick the pole (k 2 T + m 2 W − iε)/n − k for −2m χ < n − k < 0. The integral vanishes for n − k outside this range. After performing the n − k and k ⊥ integrals, we obtain the final result

B.2 Unobserved-jet function for the narrow resolution case
Here we supply the lengthy expressions for the narrow resolution jet function that were not given in [7]. The jet function of the unobserved final state in the narrow resolution regime E γ res ∼ m 2 W /m χ is defined as This is equivalent to computing the total discontinuity of the gauge boson two-point function where A B ⊥ µ is the collinear gauge-invariant collinear building block of SCET. While formally the definition appears the same as (55) for the intermediate resolution, in the present case p 2 ∼ m 2 W rather than p 2 m 2 W . The implication of this difference for the computation of the collinear integrals have been discussed in the previous subsection.
Since we are considering the χ 0 χ 0 , χ + χ − initial states and since we require a single photon in the anti-collinear final state, electric charge conservation implies that we only need to calculate the 33 component of J BC . To the one loop-order, we can write J 33 (p 2 ) as J 33 (p 2 , µ, ν) =ŝ 2 W (µ)δ(p 2 ) +ĉ 2 W (µ)δ(p 2 − m 2 Z ) + J 33 Wilson (p 2 , µ, ν) + J 33 se (p 2 , µ) . (220) where we split the result into Wilson line and a self-energy type contributions as shown in the first and second line of Figure 10, respectively . Only the Wilson line diagrams require rapidity regularization. After subtracting both dimensional and rapidity regularization poles their sum is given by where The self-energy contribution J 33 se (p 2 , µ) is expressed in terms of standard one-loop gaugeboson self-energies which can be found in [41] in the Feynman gauge. We take the fermions to be massless except for the top quark. Hence we further separate the massless fermion contribution from the massive contributions, where the second term includes the W + W − , ZH and tt loops. For the massless fermion contribution we obtain In this case we do not need to introduce star distributions, because the imaginary parts vanish below the massive thresholds indicated by the subscripts, and hence there are no singularities at 0 and m 2 Z . For convenience we collect below the explicit expressions for the gauge-boson self energies and their derivatives in the Feynman gauge. Their transverse parts are taken from [41]. The derivatives have been computed in a straightforward way.
Here m f,i is the mass of a fermion, where i indicates the generation index and f refers to the fermions within a generation. N f C is the number of fermion colors, N f C = 1 in the case of leptons and N f C = 3 in the case of quarks. The electroweak couplings are written in terms of the charge and the third SU (2) generator the previous subsection shows that the divergent terms in the integral p 2 max 0 dp 2 J 33 (p 2 ) arise from the light-fermion self-energy diagrams as these give rise to the star distributions in the one-loop result (220). Reordering this expression as we therefore focus on the last two terms. It is not necessary to perform the full Dyson resummation, resumming the Z-boson propagator insertions is enough. We obtain (dropping the argument µ from the coupling and the jet function) Note that (248) also includes the tree-level Z-boson contribution to the jet function through the last line. The terms Re Σ ZZ T (m 2 Z ) in the denominator ensure that the real part of the renormalized Z-boson self-energy vanishes at p 2 = m 2 Z as is required in the adopted on-shell scheme for the Z mass. The imaginary parts of the square brackets in (248) can be further simplified by noting that the Dyson resummation is necessary only when p 2 ≈ m 2 Z . We then obtain where denotes the tree-level decay width of the Z boson into the light fermions of the SM (all, except the top quark, masses set to zero). In deriving (249), we used the identities Re and xδ (x) = −δ(x), valid for non-singular test functions at p 2 = m 2 Z . The Dyson-resummed expression (249) can be obtained from the fixed-order expressions (220), (224) and (228) by employing the substitution rules (254)

C Soft function
In this Appendix we discuss the one-loop computation of the soft function. We start by discussing the scalar integrals for the virtual and real parts of the soft function. The final result is given by linear combinations of these integrals. We also illustrate how the rapidity divergences change between the two factorization theorems presented in the main text. Furthermore we give the inverse Laplace transforms of the resummed soft functionsŴ . The integrated soft function was defined in (59) and the index-contracted version in (60). For the calculation of the integrals and the soft coefficients we find it convenient to shift the position of the Wilson line to 0 and to perform the integration in n + y. This leads to Figure 11: Diagrammatic representation of the one-loop soft function.
Diagrammatically, the one-loop soft function is shown in Figure 11, where a single soft gauge boson attaches to any two distinct (red) dots on the external legs. In the following, we categorize the integrals according to which external legs the soft radiation attaches to. If, for example, the soft gauge boson connects the collinear (n µ − ) and the anti-collinear (n µ + ) external leg, we call this the n + n − virtual or real integral, depending on whether the soft gauge boson passes through the cut.

C.1 Virtual soft integrals
In this section we report the calculation of the relevant scalar integrals for the virtual soft correction. We define the integration measure (and implicitlyμ) as where d = 4−2 and γ E is the Euler-Mascheroni constant. For the real integrals discussed later, we also use the phase-space measure in terms of light-cone coordinates where k 2 T = −k 2 ⊥ > 0, and the delta-and theta-functions enforce n + k, n − k ≥ 0.
The n + n − virtual integral We start by analyzing the virtual n + n − integral It is convenient to proceed by first doing the contour integration in the variable k 0 . To this purpose we rewrite the integral as where If k 3 > 0 one finds four poles in the k 0 complex plane situated at ±(E k − iε), k 3 − iε and −k 3 + iε. Two of these poles are in the upper half plane and the other two are in the lower half plane. We close the integration contour in the lower half plane (notice that by doing this we pick up a factor −2πi). For k 3 < 0 we find again two poles in the upper half plane and two poles in the lower half plane. The poles in k 3 − iε and −k 3 + iε moved from the positive to the negative k 0 domain and vice versa, respectively. We obtain By summing the first terms in the two square brackets of (260) (which give the same contribution as can be easily seen by making the variable transformation k 3 → −k 3 in the second line) and after performing the k 3 and k ⊥ integrations one obtains for this part A pure imaginary part comes from the sum of the remaining terms in the upper and lower line of (260). The k 3 integral of these two terms is where the result is independent of the small imaginary part iε at O(η 0 ). After performing the k ⊥ integration and summing the two contributions we obtain the final result We perform the integration in k 0 noting that the pinched poles at k 0 = ±iε must not be picked up. These poles correspond to the potential region and are already taken into account in the one-loop contribution to the Sommerfeld effect. The integration in k is then straightforward. The result is In [6] the integrals I virt. n + n − , I virt. n + v , I virt. n − v were already computed. We find agreement except for the integral I virt.
n + n − where we find an additional term which results in the imaginary parts of (263).

C.2 Real soft integrals
The real emission contribution is extracted by applying the Cutkosky rules to the cut propagators in the previous integrals 1 We still have to keep the rapidity regulator to regulate the limit ω → 0, therefore we introduce star distributions [30], see (225) for the definition.
We first perform the integration in n − k using δ(ω − n − k), which leaves In performing the n + k integral next, the step function can be dropped as it only ensures n + k > −ω, but ω ≥ 0. Since k 2 T + m 2 W > 0, the delta-function contributes only for positive n + k, i.e. the step function does not pose a further restriction. Hence, We pulled a factor of ω η into the absolute value, as ω ≥ 0. The absolute value inside the integral forces us to consider two cases. Either m W > ω, in which case the absolute value can be dropped, as k T , m W > 0. Or ω > m W , then we split the integrand into an integral from 0 to ω 2 − m 2 W with a factor of (−1) η from the absolute value, and a second integral from ω 2 − m 2 W to ∞, where the absolute value can be dropped. To make the structure more transparent, we perform the substitution k T = k T /m W and define ω = ω/m W , which turns the previous integral into the integral over dimensionless quantities. We start with the first case ω < 1. The absolute value can be dropped and the integration results in with 2 F 1 the hypergeometric function. This is the exact result to all orders in , η. The dimensionless terms inside the curly brackets are finite in the limits ω, η → 0. Therefore the only terms involving η-poles may come from ω η−1 m −η W = δ(ω) η + 1 ω [m W ] * + O(η) in front of the bracket. The δ(ω) term in this identity requires to expand the expression in the curly brackets up to order η 1 , but allows to set ω = 0 in the function arguments. Therefore we can simplify the hypergeometric 2 F 1 function in this case. For the term involving the star distribution, we only keep the η 0 term in the bracket, which is ω independent. Therefore to order O (η, ), the result can be written as We perform the integration in n + k, n − k using the two delta-functions as for the n + n − case, and obtain Other than for the n + n − integral, the prefactor is now ω η+1 , which is finite in the limit η, ω → 0, regardless of how the limit is taken. Hence, at this point we can set η to 0. The expression is then a standard integral, that is easily solved. The result reads and is finite. To compare with the other terms, we may also replace 1 ω → 1 ω , as the integral is non-singular as ω → 0 and hence the star-distribution is equivalent to ω −1 .
The vn − real integral This integral is related to the n + n − and vn + integrals. The identity (n + · n − ) (n + k)(n − k) ensures that the integral obeys the relation I real vn − = I real n + n − − I real vn + .
The vv real integral The calculation of the vv integral follows the logic of the vn + integral. We start with the expression The two delta-functions are used for the n + k, n − k integrations. This results in the expression For the same reasons as for the vn + integral, we can expand the integrand in η and drop the O(η) terms. Hence The same result is obtained if we keep the full η-dependence and expand the hypergeometric functions that arise for the full integrals. All the real integrals except I real vv were also computed in [6], and we confirm these results.

C.3 Cut two-loop diagrams
The integrals now allow us to determine the total discontinuity of a given two-loop diagram, after summing over all cuts. The sum of all cuts is cuts = Disc(iM) = −2 Im M .
In Figure 12, we show the four possible cuts of the n + n − -diagram. For the other diagram types, we apply the same procedure. The total discontinuity for the n + n −type-diagrams is 15 Disc(iM n + n − ) = 2 Re I real n + n − − I virt.
In the n + n − -diagrams the real contribution cancels half of the virtual η rapidity divergence.
There is more than one vn − two-loop diagram. We discuss only one of them, as the others differ only via relative overall signs and prefactors: For these diagrams the rapidity divergence is completely cancelled. Following the same logic, the vn + -diagrams give The η-divergence for this two-loop diagram is the same as for the virtual integrals only, as the corresponding real integral is η finite. We observe that the left-over rapidity divergences among the two-loop diagrams are such that virt.
In the narrow resolution case E γ res ∼ m 2 W /m χ , the rapidity divergence in the sum of all virtual soft diagrams cancels the rapidity divergence of the photon jet function and the narrow resolution unobserved-jet function. The fact that the intermediate resolution case allows for real soft radiation implies that it has only half the η-divergence in the soft sector compared to the narrow resolution case. This matches precisely to the fact that the jet function for the unobserved final state at intermediate resolution and hard-collinear virtuality O(m χ m W ) does not have rapidity divergences anymore.
Using the expressions for the anomalous dimensions from the RG and RRG equations for the hard function (70), the photon jet function (76), (82), the unobserved-jet function (303), (304) and the soft coefficients of the narrow resolution case (299), (300), we can explicitly check that these consistency constraints are satisfied.

E Complete NNLO expansions
In Section 5 we showed the coefficients c (+−)(+−) defined in (138) for n ≤ 2. In this Appendix, we list the remaining coefficients relevant to the understanding of the logarithmic structure of NLL' resummation at NNLO, that is, c (n,m) IJ with 0 ≤ n ≤ 2 and 0 ≤ m ≤ 2n and I, J ∈ {(00), (+−)}. Additionally, this Appendix provides more details on the computation of the resummed cross section, as well as the fixed-order expansions. We introduce the following abbreviations (partially already given in Section 5 but repeated here for completeness): where (00)(00) = −1 − π 2 l 2 µs + l µs 2λ R l νs − 4λ R l R −