Resummed Photon Spectra for WIMP Annihilation

We construct an effective field theory (EFT) description of the hard photon spectrum for heavy WIMP annihilation. This facilitates precision predictions relevant for line searches, and allows the incorporation of non-trivial energy resolution effects. Our framework combines techniques from non-relativistic EFTs and soft-collinear effective theory (SCET), as well as its multi-scale extensions that have been recently introduced for studying jet substructure. We find a number of interesting features, including the simultaneous presence of SCET$_{\text{I}}$ and SCET$_{\text{II}}$ modes, as well as collinear-soft modes at the electroweak scale. We derive a factorization formula that enables both the resummation of the leading large Sudakov double logarithms that appear in the perturbative spectrum, and the inclusion of Sommerfeld enhancement effects. Consistency of this factorization is demonstrated to leading logarithmic order through explicit calculation. Our final result contains both the exclusive and the inclusive limits, thereby providing a unifying description of these two previously-considered approximations. We estimate the impact on experimental sensitivity, focusing for concreteness on an SU(2)$_{W}$ triplet fermion dark matter - the pure wino - where the strongest constraints are due to a search for gamma-ray lines from the Galactic Center. We find numerically significant corrections compared to previous results, thereby highlighting the importance of accounting for the photon spectrum when interpreting data from current and future indirect detection experiments.


Introduction
The discovery of the dark matter (DM) particle(s) is one of the central goals of the high energy physics program. While the Weakly Interacting Massive Particle (WIMP) paradigm with DM masses of order the electroweak scale ∼ 100 GeV has received the most attention, it is also a reasonable possibility that the WIMP could be much heavier. The canonical example is the neutral component of a new Majorana SU(2) W triplet fermion -this wino DM will be the concrete example studied here, although many of the results presented JHEP03(2018)117 below will hold for a wide class of heavy WIMPs. Assuming no other new states are present, the wino mass is the only free parameter in this model. The wino is a prototypical heavy WIMP: a calculation of the relic density for winos annihilating to electroweak gauge bosons (including the impact of the charged wino states via the Sommerfeld enhancement [1][2][3][4][5]) yields a mass of around 3 TeV. The wino as DM is motivated both from a "complete" theory perspective in the context of split supersymmetry [6][7][8][9][10][11][12][13], but it is also interesting due to its economy, i.e., minimal DM [5,[14][15][16][17].
Multi-TeV WIMPs are unobservable at the LHC: 14 TeV projected limits on winos are in the few hundred GeV range, and they will even be challenging to find at a future collider [18,19]. Furthermore, the cross section at direct detection experiments suffers an accidental cancellation between the spin-0 and spin-2 contributions, yielding a rate that is near the neutrino floor [20][21][22]. The one known channel that holds promise for detecting multi-TeV winos is via astrophysical searches for their annihilation products. Annihilation to photons could provide a very clean signal visible to ground-based air Cherenkov array telescopes [5,23,24], and constraints from the observed flux of antiproton cosmic rays can also be relevant, but require modeling of cosmic-ray propagation and backgrounds [25]. In particular, a search for line photons by the HESS experiment [26] provides a powerful constraint for thermal winos with mass near 3 TeV, although this is subject to large uncertainties from the unknown shape of the DM density profile in the inner Galaxy [23]. Furthermore, there are many upcoming experimental searches which could discover heavy WIMPs via indirect detection of gamma rays, including new data from HESS [27,28], HAWC [29][30][31], CTA [32], VERITAS [33][34][35], and MAGIC [36,37]. We would therefore like to have reliable theoretical predictions for the particle physics contribution to the cross section over a wide range of DM masses. One key feature of these ground-based experiments is that their resolution for line searches is not particularly sharp, implying that finite bin effects should be accounted for when making a precise prediction of the annihilation cross section. A main goal of the present work is to address this.
It is by now well understood that the calculation of the annihilation rate is complicated by the presence of multiple hierarchical scales, namely m W and M χ . For models with M χ m W , this separation of scales invalidates the standard perturbative expansion, introducing a number of effects that must be treated to all orders, in particular Sommerfeld enhancement, which resums terms of the form (α W M χ /m W ) k [2,3,5,38,39], and Sudakov double logarithms α W log 2 (M χ /m W ) [40][41][42][43][44][45][46]. These can be conveniently treated using effective field theory (EFT) techniques, which allow for a systematic expansion in m W /M χ 1, and the identification of universal behavior in this limit. This has attracted recent attention, resulting in calculations from different groups, with differing assumptions. Two groups [42,43,46] resummed the logarithms that appear assuming the final state was specified as γ γ or γ Z (referred to here as exclusive), while [41,44,45] calculated a resummed cross section using the operator product expansion (OPE) and assuming a γ + X final state (referred to here as inclusive). Due to these differing assumptions, distinct conclusions were reached on the importance of the logarithmically enhanced terms.
In reality, the situation is more subtle and lies somewhere in between these two extremes. Due to the finite energy resolution of the detector, the state X recoiling against the JHEP03(2018)117 detected photon, which we take to have energy E γ , is not forced to be a single electroweak boson. However, X is constrained to lie near the light cone, namely it is a jet. In this region it is well known that the standard OPE breaks down, and a more complicated factorization, describing the dynamics of the radiation within the jet, is required. Explicitly, this introduces another small parameter (1 − z) 1, where controls the distance from the endpoint, thereby further complicating the perturbative structure. In particular, large logarithms of (1 − z) appear. We will refer to these as endpoint logarithms since they become important as z → 1. The importance of these endpoint logarithms in the DM case was noticed in [45] where an attempt was made to extend the OPE based expansion beyond its region of validity into the endpoint region. 1 However, this framework did not provide a way to exponentiate these logarithms. Their resummation is one of the goals of this paper. In this paper we develop a comprehensive EFT framework to compute the photon spectrum for annihilating (or decaying) DM. We use the soft-collinear effective theory (SCET) [52][53][54], and its recent extensions developed for treating similar multi-scale problems in jet substructure, to factorize the dynamics at the scales m W (electroweak breaking scale), M χ (1 − z) (soft scale), M χ √ 1 − z (jet scale), and M χ (hard scale). In order to perform the resummation, we will need to refactorize the cross section using techniques for multi-modal field theories [55][56][57][58][59]. All large logarithms present in the cross section are then captured by renormalization group evolution between the relevant scales. The end result is a completely factorized description that allows for systematically improvable calculations of the photon spectrum. In this paper we will use this framework to compute the resummed spectrum for pure wino annihilation. The extension to Higgsinos and more general representations will be left for future work.
An example of the result from our calculation is shown in figure 1. Here we have plotted the cumulative spectrum, σ(z cut ) = 1 zcut dz dσ dz . (1.2) A value of z cut = 0 corresponds to the fully inclusive case, and z cut = 1 to the fully exclusive case. As a benchmark, we have taken the wino mass to be 3 TeV -a wider range of masses are presented below in section 6. Here we see the impact of resumming the endpoint logarithms: there is the known factor of 2.2 difference between the exclusive and inclusive calculations, and when we take z cut ∼ 0.8-0.9 (which is motivated by the HESS energy resolution), we find that the prediction falls almost half way between the inclusive and exclusive limits. This demonstrates the importance of the study presented below. The resummed cross section as a function of the experimental resolution parameter z cut for a 3 TeV wino, showing the transition between the fully inclusive (z cut = 0) and the fully exclusive (z cut = 1) cases. For z cut ∼ 0.8-0.9, as relevant for the HESS experiment, the prediction is half way between the two limiting cases, emphasizing the importance of properly treating z cut .
An outline of this paper is as follows. In section 2 we carefully review the kinematics of indirect detection, highlighting the different regions of the photon spectrum, the appropriate field theoretic techniques that are required for their description, and the differing approximations made in previous presentations. In section 3 we review the different effective theories that we will make use of in our analysis, namely non-relativistic DM effective theory (NRDM) and SCET. In section 4 we present our factorization formula for the region m W M χ (1 − z) M χ . We describe in detail the multi-step matching procedure used in its derivation, and the physical role of the different functions appearing in the factorization. In section 5 we perform the LL resummation, and derive a compact analytic expression for the resummed spectrum. In section 5.3 we show that our EFT reproduces the resummation in both the OPE region, and the exclusive endpoint by taking appropriate limits, hence tying together different results in the literature. In section 6 we present numerical results for the case of wino DM, comparing with previous results obtained using the exclusive and inclusive calculations, allowing us to demonstrate that properly accounting for the finite resolution has a numerically significant effect. In section 7 we estimate the impact of our newly derived predictions on indirect detection constraints using a simplified mock analysis of the HESS data. We conclude in section 8. Two appendices are provided: in appendix A, we provide many technical aspects of the one-loop calculations presented in the text, and appendix B demonstrates the minimal impact of photons from cascade decays (e.g. χ χ → W + W − → many γs) on our mock reanalysis of the HESS data.

JHEP03(2018)117
Guide for the reader. We anticipate that our audience's interests span from the technical aspects of the EFT-based calculation to an interest in the implications for the indirect detection experimental predictions. We therefore provide two road maps for navigating this paper, depending on the expertise of the reader. For the EFT enthusiasts, the main technical details of the factorization are presented in sections 3-5. While we have attempted to make the presentation as self contained as possible, in particular by reviewing the relevant technology, these sections necessarily assume a higher level of familiarity with EFT techniques, and are as such more mathematically intensive. These sections provide the details which yield the final prediction, but can be skipped without affecting one's big picture understanding of this work.
For the reader interested primarily in the results, and the resolution of previous differing approximations and conclusions in the literature, we recommend section 2.1 and sections 6-7. Section 2.1 emphasizes the physical differences between the different approximations previously made in the literature, and explains the necessity of pursuing our approach to derive a complete understanding for the range of parameters of interest to current and future experiments. The main results of our study are shown in graphical form in section 6, where we highlight the numerical impact of the resummation of logarithms of z cut , and compare with numerical results from previous approximations. This clearly illustrates the importance of properly including the finite resolution of the experiments. Finally, the impact of our updated numerical results on DM exclusions are given in section 7.

Kinematics for heavy WIMP annihilation
In this section, we discuss in detail the kinematics of heavy DM decay or annihilation to photons as relevant for indirect detection. We carefully analyze all relevant scales, identifying regions where large ratios of scales exist, which will give rise to logarithms that need to be resummed. This analysis will also make clear the differences between the previous studies in the literature. We will also highlight how collinear-soft modes appear in the broken theory, highlighting the distinction with the case of the naively similar B → X s γ that has been thoroughly treated in the literature (see e.g. [60][61][62][63][64]). The discussion of this section is completely independent of the details of the DM, allowing us to simultaneously consider decay and annihilation, and depends only on the kinematics of indirect detection.

Three effective field theory regimes
We consider for concreteness the annihilation of two nearly stationary DM particles of mass M χ decaying to γ + X, where the γ is assumed to be detected by the experiment. Here X denotes all final state radiation apart from the photon. The case of DM decay for a particle of mass 2M χ is identical. We use a dimensionless variable z to characterize the energy fraction of the photon E γ = M χ z , (2.1) or equivalently, Here the state X has a large invariant mass and can be integrated out. (c) The endpoint region, m X M χ . Here the measurement on the final state X constrains it to have a small invariant mass. This implies that X cannot be integrated out and must be treated as a dynamical object in the EFT. In all cases, the dashed lines dressing the annihilating DM represent the Sommerfeld enhancement.
where m X is the invariant mass of the final state X. The result of the calculation will be a differential cross section as a function of z, which will be integrated from z = z cut → 1. Depending on the value of z cut , a number of different field theoretic descriptions are required: 2 [42,43,46]: here the final state is exactly specified, either γ γ or γ Z, and we have z cut = 1. Electroweak Sudakov double logarithms, log 2 (2 M χ /m W ), appear in the perturbative expansion. See figure 2a.
• Inclusive final state ((1 − z cut ) ∼ 1) [41,44,45]: here the final state is γ + X, and the final state X is fully inclusive. This implies that m X is large, such that the state X can be integrated out using a local OPE [67]. See figure 2b.
• Endpoint region (0 < (1 − z cut ) 1): in this region, the invariant mass of the final state m X → 0 and as such it cannot be integrated out using a local OPE. The photon of interest is taken to lie along one lightcone. Then X consists of collimated high energy radiation along an orthogonal light cone, with transverse spread . The standard OPE approach is not sufficient, and a more complicated factorization theorem describing the dynamics of the soft and collinear radiation is required [68]. Deriving an analogous factorization for the case of WIMP annihilation is one of the main results of this paper. In this region, Sudakov double logarithms, log 2 (1 − z) appear in addition to electroweak Sudakov double logarithms log 2 (2 M χ /m W ). See figure 2c.

JHEP03(2018)117
We can now determine which of the above regions are most relevant to model the input photon spectrum for a search for DM lines. In principle, if the energy resolution of the detector is sufficiently precise, the appropriate cross section would only include the exclusive final state consisting of a photon and a single recoiling electroweak boson. In this case, the kinematics dictate that this condition is equivalent to requiring z 0.99 (0.9999) for M χ ∼ 500 GeV (10 TeV). The corresponding energy resolution is well beyond the capabilities of existing detectors. For example, translating the Gaussian width of the resolution quoted in the HESS line search [28] to a hard cut, would naively imply that z cut varies from 0.83 to 0.89 as M χ goes from 500 GeV to 10 TeV. This range additionally implies that we are outside the inclusive region, such that factors of log 2 (1 − z) are potentially large and resummation should be performed. We conclude then that the theory which best describes the line observations made by air Cherenkov telescopes has a state X that is recoiling against the photon with m X M χ , i.e., the endpoint region EFT. The theoretical descriptions of the matching to the exclusive region, as well as the OPE region, are also important for a complete description of the spectrum. We will see that these limits arise naturally from our endpoint EFT.

Kinematics of the endpoint region
Having determined that experimental considerations drive us to focus on the endpoint region, next we describe the relevant kinematics. This will expose the corresponding modes that will be required to construct the EFT description. These modes are shown schematically in figure 3, along with their virtualities and rapidities. Our goal in this section is twofold. First, this discussion will motivate the EFTs introduced in section 3. Second, it will allow us to provide context and highlight the new features of the factorization needed here in a physical manner, motivating the technical discussion of section 4. The later sections will then provide a comprehensive mathematical treatment, to complement the simple picture that follows from kinematic arguments.
We begin with the kinematics of the initial state, namely the annihilating DM. The DM in the halo has a typical velocity v ∼ 10 −3 , so a non-relativistic description is appropriate. The DM will be modeled as heavy sources (in analogy with heavy quark EFT or nonrelativistic QCD) emitting ultra-soft radiation, as shown in figure 3a. There is one well known complication in the heavy mass limit. Winos carry electroweak charge such that the Sommerfeld enhancement due to the exchange of electroweak gauge bosons must be included. This can be appropriately accounted for in the non-relativistic DM (NRDM) EFT by including the relevant potentials, see section 3.1. A feature of the NRDM EFT is that it allows a factorization of the Sudakov corrections from the Sommerfeld effects.
The final state is more complicated, and a full characterization will require a multimodal EFT. Recapping the discussion above, as z cut → 1 the final state consists of both a jet of collimated energetic particles and wide angle low energy radiation. As is well known, this can be captured by SCET. However, due to the multi-scale nature of the problem, we will show that additional modes, illustrated in figure 3, will be required to fully factorize all the logarithms. The origin of the multi-modal structure, and its complexity compared to that seen in previous approaches to heavy WIMP annihilation, can be understood from JHEP03(2018)117  kinematic arguments. Specifically, logarithms appear due to two types of phase space restrictions: • Kinematic restrictions on final states of massless particles: 3 These include kinematic restrictions via event shape observables, such as thrust, or restrictions from kinematics that force one into an endpoint region, as in B → X s γ, and have been discussed above. EFT descriptions in these cases typically involve three scales: the hard scale, which in our case will be M χ ; the scale of the transverse momenta of particles in the jet (whose modes are called collinear), namely M χ √ 1 − z; and the energy scale of soft radiation, namely M χ (1 − z). This class of problems is well understood and can be treated using SCET I , discussed in section 3.2. The radiation in the final state is factorized into energetic modes, referred to as collinear (c), which comprise the dynamics of the jet, and wide angle low energy radiation, referred to as ultrasoft (us). Decomposed into light cone coordinates (n · p,n · p, p ⊥ ) (see eq. (3.7)), along the direction of the jet, these modes have momentum scaling as 4 Here we mean massless in perturbation theory, as relevant for scales appearing in logarithms in the weak coupling expansion. Other mass scales can appear non-perturbatively, for example, hadron mass effects in QCD event shapes have been studied in [69,70]. 4 Note that here and throughout the text, when we describe the scaling of modes we indicate only the parametric scaling as a function of the relevant scales in the problem, namely Mχ, mW , and 1 − z. Any O(1) numerical factors do not modify this scaling, and are therefore neglected.

JHEP03(2018)117
• Exclusive final states of massive particles: these include the classic massive Sudakov form factor [71], and more recently, exclusive electroweak production [72][73][74][75][76][77], and the exclusive approximation for DM annihilation discussed above [42,43,46]. Here there are two relevant mass scales, namely the hard (h) scale M χ , and the scale of the massive boson, m W . Problems of this type can be treated using an SCET II theory, discussed in section 3.2. The relevant modes in the effective theory are collinear (c) and soft (s) modes. Decomposed into light cone coordinates (see eq. (3.7)), along the direction of the jet, these modes have momentum scaling as Note the distinction in scaling between the ultrasoft and soft modes. While in this case the collinear and soft modes are at the same virtuality p 2 ∼ M 2 χ λ 2 , they are separated in rapidity. 5 This explains the appearance of the rapidity axis in figure 3b.
The annihilation of WIMP DM in the endpoint region is a more complicated problem, since it involves the physics of both types of restrictions. There is both a constraint on the final state radiation, as well as the presence of the mass scale of the electroweak bosons and the measurement of just the photon state from among the SU(2) × U(1) gauge bosons. Indeed, we will find that all the scales (in both rapidity and virtuality) present in both individual cases will appear. This is illustrated in figure 3b, which shows the modes that live at each of these mass and rapidity scales. We will show how to factorize the dynamics at each of these scales when large hierarchies are present, thereby facilitating resummation. The final form involves a component where the gauge boson can be treated as massless, so that the scale is set by the final state kinematic restriction, and a component where the relevant scale is m W . For example, the description of the final state jet will be split into a massless jet function, described using standard techniques in SCET I , as well as a function describing the dynamics at the scale m W , using SCET II .
In addition to these SCET I and SCET II ingredients, we will show that an extra mode is required to achieve the fully factorized result. This mode has a virtuality µ 2 ∼ m 2 W , but it has a large momentum component along the direction of the recoiling photon of size M χ (1 − z) (the momentum scale of the soft function): .
In the case that both M χ (1 − z) M χ and m W /(M χ (1 − z)) 1, these modes are neither (ultra)soft, or collinear, i.e., they do not appear in either SCET I or SCET II EFTs, but are instead an example of collinear-soft modes, see section 3.2. Our factorization formula allows for the separate treatment of these collinear-soft modes, which allows us to resum all large logarithms, but also ensures continuity of the cross section as we move away from the endpoint region, where these modes are no longer distinguishable from the standard 5 We will typically use a dimensionful rapidity, ν, as in figure 3b. This should be thought of in analogy with the dimensional regularization scale, µ, and is introduced in section 3.2 where we discuss the regularization of rapidity singularities.

JHEP03(2018)117
soft modes. It is the simultaneous presence of the scales M χ (1 − z) and m W that gives rise to the presence of these collinear-soft modes -they would not appear if only a subset of the scales were present. 6 The structure of the results presented below shares similarities with the factorization formulae for jet substructure observables, where a measurement in addition to the mass has been performed [55][56][57][58][59][78][79][80].
The complete description of the final state therefore combines the SCET I collinear and ultrasoft modes with the SCET II soft and collinear modes in the direction of the jet, along with the collinear-soft modes describing additional radiation along the direction of the photon. Each of these will yield distinct functions in our factorization formula eq. (4.1), implying that each of these functions has a clear physical origin in terms of the scales of the problem. This seemingly complicated description is in fact a significant simplification, since the description of the dynamics at any one of these scales has been reduced to its elemental form. In the next section, we will introduce the EFT ingredients, and in section 4 we give the technical details of the factorization.

Review of relevant effective field theories
In this section we briefly review the different EFTs that we will use, primarily to establish our notation. Our use of non-relativistic (NR) field theories will be standard in the context of QCD [81][82][83] (for reviews, see [84][85][86]), and will focus on aspects relevant for annihilating DM (for applications of NRDM EFT to the scattering of DM with nucleon targets, see [20,21,87,88]). As we review SCET, we will highlight necessary extensions that are perhaps less familiar.

Non-relativistic dark matter effective theory
In the NRDM EFT, large fluctuations of the DM field χ about a particular velocity v are integrated out. The non-relativistic DM is described by a field χ v with a label velocity v, just as in heavy quark EFT [89,90]. Here v is a dimensionless four vector describing the velocity of the DM, which for concreteness we will take to be v = (1, 0, 0, 0). The freedom in the choice of v is represented in the EFT as a symmetry known as reparametrization invariance [91,92]. The dynamics of χ v describe the residual fluctuations of the heavy state, as in non-relativistic QCD. The EFT captures the interactions of the non-relativistic particles whose momenta p µ = (E, p ) scale as soft (M χ v, M χ v), ultrasoft (M χ v 2 , M χ v 2 ), and potential (M χ v 2 , M χ v). The ultrasoft modes describe radiation, while the soft modes give rise to the running of potentials.
The leading power interactions of the heavy DM particle(s) with the ultrasoft radiation can be eliminated using a field redefinition χ Here we have argued for the existence of collinear-soft modes based only on kinematics. The fact that these modes are actually required is also related to the fact that there are external states with electroweak charges, as will be discussed in section 4.2.

JHEP03(2018)117
where P denotes path ordering, g is the relevant gauge coupling, and T a (r) is the generator for the DM representation r. Furthermore, soft radiation is not required at the order to which we work. This implies that all dynamical radiation in NRDM is completely captured by Wilson lines along the directions of the heavy particles, greatly simplifying the field theory treatment.
After decoupling the soft radiation, the leading power Lagrangian is given by which describes the interactions of the heavy particles as the sum of a kinetic and potential term. The potentialV describes potential exchanges of the W, Z, γ, and its explicit form can be found in ref. [4]. Note that going to higher orders and powers is well understood in the context of NRQCD (see e.g. refs. [93,94]). The dynamics of the heavy particles are governed by low energy matrix elements evaluated with the above Lagrangian. Since this is a non-relativistic description, the number of heavy particles is fixed, and there exists an associated Schrödinger equation. These low energy matrix elements give rise to the Sommerfeld enhancement, which must be included when computing the DM cross section. We will therefore briefly review the structure of the low energy matrix elements and the Sommerfeld factors.

Sommerfeld factors
Since we have chosen to work with pure wino DM, the model includes a Majorana fermion DM candidate χ 0 , and an electrically charged fermion χ ± . For the calculation of the Sommerfeld factors, we include a mass splitting, that is neglected when performing the Sudakov resummation. Including this splitting is important as it plays a role in determining the positions of the Sommerfeld resonances. For winos, electroweak corrections yield a mass splitting δ ≡ M χ ± − M χ 0 164.4 MeV [95]. In our formalism, the Sommerfeld enhancement will be captured by low energy matrix elements of the heavy annihilating particles. As discussed in section 4 where we derive the factorization formula, the following matrix elements appear where T denotes transpose, σ 2 is the second Pauli matrix, and the external state is given by the S-wave combination (χ 0 χ 0 ) S . Here the color indices a, b, a , b = 1, 2, 3, and we have the usual relations χ 0 = χ 3 and χ ± = (χ 1 ∓ iχ 2 )/ √ 2. In terms of the charge eigenstates, we will find that the relevant components of F a b ab are where the Sommerfeld enhancement is captured by the factors s 00 and s 0± , which must be evaluated non-perturbatively. In practice we do this by numerically solving the associated Schrödinger equation. We summarize some of the most important aspects here; a

JHEP03(2018)117
detailed discussion can be found in appendix A of [23]. For other detailed studies of both phenomenological and formal aspects of Sommerfeld enhancement, we refer the reader to refs. [96][97][98][99][100]. The first step in solving for the Sommerfeld factors is to compute a wavefunction ψ i j , where the index i labels the asymptotic state and j is the component index for the resulting solution, and the indices i, j = 1, 2 refer to the (00), (+−) states respectively. A discussion of the relevant boundary conditions can be found in ref. [23]. Once the solutions ψ have been obtained, the Sommerfeld enhancement matrix is (3.5) In practice, one must choose a velocity when computing s ij . As is well known, the Sommerfeld enhancement saturates at low velocities, and we have checked that this occurs for the range relevant for DM annihilations, i.e., v 10 −3 , for the wino mass range of interest. Therefore, we can neglect any velocity profile dependence, and treat all velocity dependence as constant for the parameter range of interest.
Once we know s ij , using eq. (3.4) we can then determine the relevant components of F a b ab given in eq. (3.3). From this point, the annihilation cross section can be computed as whereσ a b ab (z cut ) denotes the resummed perturbative cross section as a function of z cut , whose computation is the subject of this paper (see eq. (4.1) below). As a final comment, we note that we have glossed over the fact that we will be working in a theory with a spontaneously broken gauge symmetry, as opposed to standard NRQCD. There will several manifestations of this fact. First, and most trivially, it impacts the Sommerfeld enhancement calculation, as well as the color algebra, due to the identification of a color index for the external photon. More non-trivially, a significant portion of this paper (see in particular section 4) will relate to the refactorization of the function describing wide angle soft radiation, including that from the incoming DM particles. This is required, since m W introduces another scale for the soft radiation in addition to that imposed by the final state measurement.

Soft-collinear effective theory
Soft-Collinear Effective Theory (SCET) [52][53][54] will provide the framework for describing radiation in the final state. SCET describes the dynamics of soft and collinear radiation in the presence of a hard scattering. While originally developed for applications to QCD with massless gauge bosons, the formalism was extended to the electroweak sector with massive gauge bosons in [72][73][74]. In what follows, we will provide a brief review of the features of SCET that will be used for our heavy DM annihilation process (along with a few more general comments).

Modes, fields, and Wilson lines
SCET is a theory of both soft and collinear particles. Collinear particles have a large momentum along a particular light-like direction, while soft particles have a small momentum, and no preferred direction. For each relevant light-like direction, we define two reference vectors n µ andn µ such that n 2 =n 2 = 0 and n·n = 2. The typical choice of n µ = (1, 0, 0, 1) andn µ = (1, 0, 0, −1) will be used below. The freedom in the choice of n, as in the case of v for non-relativistic EFTs, is represented in the EFT through a reparameterization invariance [101,102]. Any four-momentum p can be decomposed with respect to n µ as The SCET expansion is defined by a formal power counting parameter λ 1, which is determined by the measurements or kinematic restrictions imposed on the radiation. Then the momenta for the different particles in the EFT scale as where Q is a typical scale of the hard interaction. A theory with collinear and ultrasoft modes is typically referred to as SCET I , while that with collinear and soft modes is referred to as SCET II [103]. 7 In order to expand the full theory fields around a particular direction, the momenta are decomposed into labelp µ and residual k µ components Then for a collinear particle,n ·p ∼ Q andp ⊥ ∼ λQ, while k ∼ λ 2 Q describes small fluctuations about the label momentum. EFT modes with momenta of definite scaling are obtained by performing a multipole expansion of the full theory fields. SCET involves independent gauge bosons 8 for each collinear direction A n,p (x), which are labeled by their collinear direction n and their large label momentump, as well as (ultra)soft gauge boson fields A (u)s (x). Independent gauge symmetries are enforced for each set of fields. Overlap between different regions is removed by the zero-bin procedure [106]. This ensures that there is no double counting of momentum regions. The leading power SCET Lagrangian takes the form Here L (0) hard contains the hard scattering operators and is determined by an explicit matching calculation. The Lagrangian L (0) describes the universal leading power dynamics of the 7 In the presence of Glauber modes, soft modes are always required to run the Glauber potentials [104,105]. Whether or not ultrasoft modes are required depends on the physical observable in question. 8 The standard formalism also incorporates collinear scalars and fermions as well. These are not required for the calculation presented here, so we will not discuss them.

JHEP03(2018)117
soft and collinear modes and can be found in refs. [52][53][54]. Finally, L G is the leading power Glauber Lagrangian [104], which describes the leading power coupling of soft and collinear degrees of freedom through potential operators. We will not need to consider it in this paper.
Hard scattering operators involving collinear fields are constructed out of products of collinear gauge invariant fields [52,53]. The gauge invariant gauge boson operator is given by Here D n⊥ is the collinear gauge covariant derivative, and W n is a collinear Wilson line 9 where P µ is an operator that returns the label momentum. The collinear Wilson line, W n (x), is localized with respect to the residual position x so that B µ n⊥ (x) can be treated as local gauge boson fields from the perspective of the ultrasoft degrees of freedom. For the leading power calculation presented here, ultrasoft and soft fields will not appear explicitly in our hard scattering operators, other than through Wilson lines via the field redefinition which is performed in each collinear sector. For a general representation, r, the ultrasoft Wilson line is defined by 10 where as before P denotes path ordering. This so-called BPS field redefinition has the effect of decoupling ultrasoft and collinear degrees of freedom at leading power [111]. We will also need soft Wilson lines, Note that when the label momentum is large compared to the virtuality of the EFT modes, it is convenient to use a mixed position/momentum space representation space Wilson line, where the label is in momentum space and the residual fluctuations are in position space. Otherwise, Wilson lines will be written in position space, e.g. eq. (3.1). It is also possible to formulate SCET entirely in position space, see e.g. refs. [107,108], although we will not use the position space formalism here. 10 Here we give the explicit result for an incoming Wilson line. Depending on whether particles are incoming our outgoing, different Wilson lines must be used. When done correctly, the BPS field redefinition accounts for the full path of the particles [109,110].

JHEP03(2018)117
Finally, the refactorization of the soft sector (see section 4.2.3 below) will require the inclusion of collinear-soft modes from SCET + [55][56][57][58][59]. Collinear-soft modes have both a collinear and soft scaling where λ andλ are distinct power counting parameters. Such modes first appeared in calculations of jet substructure when multiple simultaneous measurements are made on a jet [55][56][57][58][59]. This introduces additional scales, implying the need for both λ andλ. For contrast, the measurement of a single observable, such as the mass of a jet, only introduces a single scale; the mass can either fix the angular spread of the mode, resulting in a collinear mode, or it can fix the energy of the mode, resulting in soft or ultrasoft modes, but it cannot fix both, as required for collinear-soft modes. In our case, the collinear-soft modes will arise due to the presence of both the mass scale of the final state m X , and the mass scale of electroweak symmetry breaking m W . Our study provides a new application of collinear-soft modes.
Since the collinear-soft modes arise from a refactorization of the soft sector, they couple eikonally and their interactions can be absorbed using additional Wilson lines defined as and This notation is chosen to reflect that the X Wilson lines will arise from a BPS field redefinition, similar to the Y Wilson lines in SCET I (and X precedes Y in the alphabet), and the V Wilson lines are generated by integrating out interactions with particles in then direction, similar to the W Wilson lines that accompany the collinear fields (and V precedes W in the alphabet). As with (ultra) soft fields, at the order to which we work, collinearsoft fields will appear only in Wilson lines. For example, they will arise from the BPS field redefinition, which allows the all orders decoupling of interactions between collinear-soft and collinear particles. This is identical to the transformation in eq. (3.13) but with a collinear-soft Wilson line. For a more detailed discussion of the BPS field redefinition for collinear-soft fields, see [55].

Renormalization group evolution
SCET allows for the resummation of large logarithms through the renormalization group (RG) evolution of matrix elements of collinear, (ultra)soft, collinear-soft fields. Since we will use both SCET I and SCET II , this RG evolution can be either in virtuality, µ, or JHEP03(2018)117 rapidity, ν [112][113][114]. We use the regulator of [113,114], modifying the Wilson lines as Here ν is a rapidity scale, analogous to µ in dimensional regularization, η is the regulating parameter, and P z returns the z-component of the label momentum. This allows us to define a dimensional regularization-like RG in terms of ν.
Here ω is a formal bookkeeping parameter which satisfies For convenience, we set ω = 1 throughout our calculations since it can be trivially restored. Rapidity divergences for the collinear-soft modes will also be regulated with the appropriately modified versions of eqs. (3.19) and (3.20).
In our factorization, we will encounter functions that satisfy both multiplicative and convolutional renormalization group equations. For a function F (µ, ν) which is renormalized by a multiplicative factor Z F (µ, ν), we have Convolutional renormalization in a variable τ takes the form giving rise to the RG equations where the anomalous dimensions are given by

JHEP03(2018)117
Convolutional RG equations are most easily treated in a conjugate space (we will use Laplace space below), in which they are multiplicative. The RG evolution can be used to run functions from their natural scale, where all large logarithms are minimized, to an arbitrary scale. The independence of the RG path is guaranteed by the fact that the anomalous dimensions sum to zero, schematically where the sum is over the functions F that appear in the factorization formula, along with the fact that evolution in µ and ν commutes: The consistency of the anomalous dimensions will provide a strong check on our calculation. We will use the path independence to choose a particularly simple path to resum all large logarithms in the EFT, see figure 6 below.

Factorization formula for the endpoint region
In this section, we present the factorization formula for the endpoint region of heavy WIMP annihilation -this is one of the main results of this paper. We focus here on the shortdistance component of the cross section, denotedσ(z cut ) in eq. (3.6). As discussed below, the long-distance contributions, i.e., the Sommerfeld enhancement, also arise naturally from the factorization of the matrix elements presented in this section; we refer the reader to section 3.1.1 for the details of how these factors are (numerically) computed. In section 4.1, we present the factorization formula, and discuss each of its components in turn. This section is aimed at readers without a technical EFT background, and as such emphasizes the physical content of each ingredient. In section 4.2, we provide the technical discussion of the multi-stage matching used to derive the factorization formula, emphasizing the operator definitions for the functions and key aspects of the refactorization. Tree level and one-loop results for all functions in both the intermediate and final EFT, as well as details of the calculations can be found in appendix A.

Factorization overview
The main result of this section is a factorization formula for the photon spectrum in the endpoint region. We find that the differential cross section for the heavy WIMP annihilation χ χ → γ + X factorizes in the limit that z → 1 as where z is defined in eq. (2.1), and we use ⊗ to denote a convolution between the functions in the second line, as explained in detail below. Hereσ denotes the short-distance component

JHEP03(2018)117
Increasing Virtuality H H Figure 4. A schematic of the multistage matching procedure used to derive the factorization formula for heavy WIMP annihilation in the endpoint region. The jet and soft functions appearing in the first stage of matching are refactorized into components that depend either on m W , or on the phase space restriction implemented by z.
of the cross section in eq. (3.6) with suppressed initial/final state indices. The indices are to be contracted with the matrix element F a b ab in eq. (3.3). This function also arises naturally when considering the factorization of the cross section, but to keep our discussion focused on the Sudakov factors, we will not consider F a b ab in this section. When we present the final cross section results in section 5.2, F a b ab will be included. The LL superscript indicates that this factorization as written is only true for the leading logarithmic contributions. Beyond this order additional functions are required, as will be described in this section. The iterative matching procedure used to derive this result is shown schematically in figure 4. In the first stage, we match onto a standard SCET theory, leading to the standard factorization into functions that describe the underlying hard scattering (H), the collinear radiation along the jet (J n ) and photon (J γ ) directions, and soft radiation (S ). In the second stage, we match onto a (electroweak symmetry breaking) theory with massive soft and collinear modes. In particular, this manifests as a refactorization of the soft function S into the functions H S , S and C S , and of the jet function J n into the functions H Jn and Jn -these additional functions are described below.
The final EFT description consists of a collection of independent sectors, each corresponding to the functions appearing in the factorization formula eq. (4.1). The procedure for factorizing the full cross section into these functions is illustrated in figure 4. The interpretation of each of the functions is discussed in the following, which is organized by the characteristic scale µ for these sectors. In particular, we separate it into two classes of functions, namely those that depend on m W , and those that do not.

JHEP03(2018)117
The first class of functions depend on scales far above the electroweak scale, µ m W , and are thus independent of electroweak symmetry breaking effects.
• H Jn (M χ , 1 − z, µ) describes collinear radiation along the jet direction with virtuality µ ∼ M χ √ 1 − z such that it contributes to the final state mass.
such that it contributes to the final state mass.
The second class of functions encode electroweak symmetry breaking effects, and have µ ∼ m W , so that the gauge fields are treated as massive. Additionally, these functions all depend on a rapidity renormalization scale ν.
• J γ (m W , µ, ν) describes the final state photon, and results purely from modes with energy E γ and virtuality µ ∼ m W . This function receives only virtual corrections, since the final state is exactly specified.
• S(m W , µ, ν) describes homogenous soft radiation with virtuality µ ∼ m W such that it does not contribute to the final state mass.
describes radiation that is simultaneously soft and collinear to the photon direction. The momentum for this radiation has collinear scaling, virtuality µ ∼ m W , and contributes to the final state mass.
• Jn(m W , µ, ν) describes collinear radiation along the jet direction with virtuality µ ∼ m W such that it does not contribute to the final state mass.
This full factorization simultaneously involves functions from NRDM, SCET I , SCET II , and SCET + , and resummation requires RG evolution in both virtuality and rapidity. For the analysis here, we will be interested in resumming only the leading logs (LL). Our approach to the factorization persists at higher logarithmic order. However, as written, the refactorization of the soft function S is only valid at LL order. The origin of this effect, as well as the mechanism for disentangling these scales, is akin to the case of non-global logarithms (NGLs), and is discussed in section 4.2.3.
While we will present the factorization formula using the concrete example of an SU(2) W triplet of Majorana fermions, this choice merely affects the particular spin and charge structure of the operators involved, and as such the main features of the factorization and the relevant modes are universal. The same factorization will also apply, e.g. to the annihilation of heavy SU(2) W doublets or the decay of a heavy dark bound state [115]. Furthermore, some of the structure is generic to situations where event shape observables are measured on jets of massive radiation, and thus variants of eq. (4.1) may find applications for future high energy colliders [116,117].

Multi-stage matching
In this section, we discuss the derivation of the factorization formula given in eq. (4.1). In section 4.2.1 we present the first stage of matching, including the structure of the hard scattering operators, the factorization of the Hilbert space and measurement function for soft and collinear modes, and the matrix element definitions of the functions. In section 4.2.2 and section 4.2.3 we present the details for the second stage of matching, namely the refactorization of the collinear and soft sectors. For the soft sector, we give a detailed discussion of the relevant soft and colinear-soft modes.

Soft-collinear factorization
We begin by determining the hard scattering Lagrangian in SCET, denoted by L hard in eq. (3.10). This is done through matching the full theory consisting of the Standard Model and an SU(2) W triplet of Majorana fermions onto SCET, and is identical to the fully exclusive case [42,43,46]. The Lagrangian describing the hard scattering is with the Wilson line structures obtained through the BPS field redefinition. The Wilson coefficients C r are IR finite, and independent of the scale m W . Performing a tree-level matching at the scale µ ∼ M χ , we find The C r (µ) encode the underlying hard scattering process and determine the hard function H(M χ , µ) appearing in our factorization formula, as will be defined in eq. (4.9). Together with L dyn in eq. (3.10), the hard scattering operators in eq. (4.2) describe the annihilation at scales µ M χ . The factorization formula for the cross section for χ χ → γ +X depends on the squared matrix elements of these hard scattering operators. For contrast, in the exclusive case there are only virtual contributions, and thus the factorization can be done at the level of the amplitude [42,43,46]. In the present analysis, there are both real and virtual contributions that are sensitive to m W as well as the scales imposed by the endpoint restrictions though z. These low-energy dynamics are not yet factorized at this stage.
First, we consider the factorization of the Hilbert space for the final state |X . Since the soft and collinear modes are decoupled, the final state can be written as (4.5)

JHEP03(2018)117
Next, we expand out the contributions to the final state mass m 2 X , These measurement operators can be written in terms of the energy momentum tensor of either the full or effective theories [118][119][120]. Here their role will simply be to return the value of the observable for a particular perturbative state in momentum space.
With the above ingredients, we can algebraically manipulate the cross section into a factorized form involving matrix elements of either soft or collinear fields. These matrix elements will be coupled together both through color indices and the convolutions that are present as a result of enforcing the measurements. This procedure is standard (see, e.g. the review [121]) and we simply give the final result. At the first stage of matching, the differential cross section with factorized dynamics in SCET is given in terms of the hard function H, the jet functions J n and J γ for X and the photon respectively, and the soft function S as where we have suppressed the color indices and the dependence on the RG scales µ and ν for simplicity. As in eq. (4.1), we have used ⊗ to denote the convolution in z. The convolution arises due to the fact that the total invariant mass of the final state is a sum over the soft and collinear sectors, see eq. (4.6).
The functions labeled with a superscript prime are those that require further factorization. Note that the J n and S functions still depend on both the m W and (1 − z) scales. This complication did not occur for the fully exclusive case, where the above factorization JHEP03 (2018)117 was sufficient since there is no intermediate scale (1 − z). The refactorization of the jet and soft functions will be discussed in section 4.2.2 and section 4.2.3.
Next, we provide field-theoretic definitions for the functions appearing in eq. (4.8). The hard function is defined in terms of the Wilson coefficients of the hard scattering operators in eq. (4.2) as The soft function is a vacuum matrix element of the soft Wilson lines Y r in eq. (4.2), where the color indices are suppressed, T andT denote time ordering and anti-time ordering respectively, and the Y r factors are the products of Wilson lines defined in eq. (4.3). The components of the soft function with explicit color indices are where the color indices are explicit, but we have dropped the arguments and scale dependence of the functions for simplicity. Here, as well as in the expressions below, we keep the time ordering convention and the dependence on x = 0 implicit. Note that the color index 3 corresponds to the photon final state. The indices i, j in the hard and soft functions span the space of the operators given in eq. (4.2) and are contracted with each other as H ij S ij . To reduce the number of indices appearing in later formulas, we introduce the following notation: The jet functions for the recoiling jet X and the photon are color-singlet matrix elements of collinear fields. Explicitly, we have where P ⊥ returns the perpendicular component of the label momentum. As discussed above, this is the final form for J γ , but the jet function for X will require further factorization -we turn to this in the next section.

Refactorization of the jet sector
As currently formulated, the jet function J n in eq. (4.8) results from dynamics at both the scale M χ √ 1 − z and the scale m W . To be able to resum logarithms of m W /(M χ √ 1 − z ), which can become large as we move towards the endpoint, we must factorize these two scales. This factorization is similar to that performed in the fully exclusive case, where one is separating M χ from m W using a hard matching coefficient that is independent of the IR scale m W , along with jet and soft functions which describe the dynamics at the scale m W . Here we will write the jet function J n (M χ , 1 − z, m W , µ, ν) as a hard matching coefficient H Jn (M χ , 1 − z, µ), and a jet function Jn(m W , µ, ν).
The collinear state, X c , factorizes into two types of collinear modes as where c z is in the Hilbert space containing the collinear modes that are sensitive to the measurement enforced as a function of z, while c W is in the Hilbert space that contains the modes with mass m W . This follows from the same logic as the standard hard-collinear factorization. Here the c z modes which contribute to the jet mass measurement have the standard scaling for an SCET I collinear mode associated with the mass measurement, The modes sensitive to the m W scale are standard SCET II collinear modes at the scale m W , with scaling and do not contribute to the mass of the final state at leading power. The factorization of the measurement function is trivial since, at leading power, the low-energy collinear modes have an invariant mass p 2 , and do not contribute to the mass of the final state. We therefore only have The separation of collinear modes through eqs. (4.14) and (4.17) allows us to fully factorize the jet function as This factorization is a power expansion in m W /(M χ √ 1 − z ). The matching coefficient H Jn can be evaluated in the unbroken theory with massless electroweak bosons, and is IR finite due to the mass measurement. The dependence on the electroweak scale is completely captured by the function Jn(m W , µ, ν).

Refactorization of the soft sector
Next, we turn to the refactorization of the soft function S . The goal is to have separate EFTs for the dynamics at scales µ ∼ M χ (1 − z) and µ ∼ m W . Comparing to the discussion of the jet refactorization in the previous section, the physics of the soft sector is more interesting, as logarithms due to collinear-soft modes appear.
Consider the possible classes of soft modes with virtuality µ 2 ∼ m 2 W . The virtuality of the soft modes with scaling p S ∼ M χ (1 − z)(1, 1, 1) can be lowered uniformly to yield modes with p S ∼ (m W , m W , m W ). When acting on these states, the measurement function in eq. (4.10) can be expanded as We conclude that these soft modes do not contribute to the measurement, which allows a simplification of the operator structure. As an explicit example, the soft functions S 1 and S 2 become where we have used the unitarity of the Wilson lines. These new functions S i are now independent of m W . Physically, the simplification (collapse) of the Wilson lines occurs because the measurement operator has been expanded away, implying that the refactorized soft functions are now inclusive. However, we are still specifying the photon as the final state, and therefore violate the assumptions of the Bloch-Nordsieck [122] or KLN [123,124] theorems, as originally pointed out in [125][126][127]. This explains why the Wilson lines in S 2 do not completely simplify, as compared to S 1 where the Wilson line dependence has collapsed to the unit operator leaving behind only color and kinematic factors. It is clear from the collapse of the Wilson lines that the modes p S are not sufficient to complete the picture. In particular, the divergences associated with m W , for example in S 1 , should be reproduced after factorization, but the function S 1 in eq. (4.20) does not have such a divergence. Interestingly, however, there is a second possibility for lowering the virtuality of the soft modes down to the scale m W : keep their momentum component along the photon direction fixed, but decrease their angle (increase their collinearity) with respect to the photon. These modes are shown schematically in figure 5. Such modes then have the scaling . modes discussed in section 3.2, which arise from the simultaneous presence of the two scales M χ (1 − z) and m W . These arguments imply that the Hilbert space of the soft sector factorizes into soft modes with uniform scaling and collinear-soft modes as (4.23) The soft modes do not contribute to the measurement, while the collinear-soft modes are sensitive to a measurement function The most interesting aspect of these collinear-soft modes is that they contribute to the measurement of the final state mass through their large component, which is independent of their virtuality. To our knowledge, this type of collinear-soft mode has not previously appeared in the literature. For example, in the case of thrust [128] or other SCET I event shapes, the definition of the measurement guarantees that it is always the small component of the momentum of a particle that is measured.
Using the measurement function in eq. (4.24), the Wilson lines that make up the collinear-soft function do not collapse, but are instead expanded assuming the momentum scaling for the collinear-soft modes. Since the collinear-soft modes are boosted along the photon's direction n, the v andn Wilson lines appear to collapse down to then direction. The collinear-soft function is therefore given as a product of Wilson lines where the X and V Wilson lines were defined in section 3.2, and implicitly include rapidity regulators. We have suppressed color indices for simplicity. Explicit expressions with color JHEP03(2018)117 indices will be given below. To regulate rapidity singularities for the collinear-soft Wilson lines, we do not expand the regulator, using the full |2 k z | −η dependence. Performing the naive power expansion of the regulator yields unregulated rapidity divergences in the collinear-soft sector. This choice of regulator defines the zero-bin structure [106] of the collinear-soft sector, and we find that non-trivial zero-bins are present, which must be correctly incorporated to remove overlap. This is described in more detail in appendix A. Strict power counting can be preserved by introducing a boost parameter β, and using the regulator |2 k z | −η → |k + + β k − | −η [104].
Having discussed the modes that are required to describe the physics at the scale m W , we next explain how to refactorize the soft function into a matching coefficient that describes the dynamics at the scale M χ (1 − z), and a soft and jet function that describe the dynamics at the scale m W . This is more complicated than for the jet function. The complication emerges due to the existence of a hierarchy in energy but not in angle for the homogeneous soft modes that live at the scales m W and M χ (1 − z). Hence, any emission at the scale M χ (1 − z), which can be at an arbitrary angle, eikonalizes from the perspective of the emissions at the scale m W , and is described as a new Wilson line source. In this way, an infinite number of operators is generated in the matching (although only a finite number appear at any order in α W ). This situation is familiar from the case of NGLs [129], where there exist multiple hierarchical soft scales. Due to the generation of these new sources, the resummation of NGLs is governed by the non-linear BMS equation [130]. In the present case, however, the measurement function for the modes at the scale m W is expanded, and what is generated are Bloch-Nordsieck or KLN violating NGLs. We are not aware of these appearing previously in the literature. While it is possible that these take a simple form, or completely cancel, they first contribute at NLL order. Here we restrict ourselves to LL accuracy, and so we will not discuss this higher order structure any further. We leave the study of them using existing formulations of NGLs in factorization [58,[131][132][133][134] for future work.
At LL order, we do not need to consider the generation of additional Wilson lines in the matching. Nevertheless, the general structure of the refactorized function can become complicated since four Wilson lines appear in each of the soft and collinear-soft functions, and mixing between these color structures can be generated beyond tree-level. In the most general case, the refactorization takes the form .
This refactorization, along with the scales of each of the functions, is shown in figure 5. The functions C S and S each carry eight color (triplet) indices. Two of these sixteen color indices are identified as carrying the quantum number of the photon, and the rest are contracted as to leave the overall indices aba b , which are contracted with the initial state wavefunction factors. In eq. (4.26), we are using the notation introduced in eq. (4.12); the index i enumerates the color structures in the soft function before refactorization, i.e., i = 1, 2, 3. The index j sums over the color structures in the soft function after refactorization.

JHEP03(2018)117
Instead of writing down a complete basis, we construct the color structures explicitly from the top down by explicitly refactorizing the soft function S . This requires us to supplement the operators written in eq. (4.11) above with those that appear at one-loop, to ensure that the RG closes. Fortunately, only a limited basis of color structures is required at this order. The color structures are derived in appendix A.2. Here we simply state the results for the refactorization of the soft functions. We denote the combined collinear-soft and soft functions asS and they arẽ (4.28) Here we have made the color structure explicit, but we have dropped the arguments and scale dependence of the functions for simplicity. The collinear-soft function reproduces the m W dependent IR divergences of the soft function. Additionally, for the RG to close we will need the following operator which has a vanishing tree-level matching coefficient, but will appear in the mixing that results as we RG evolve the functions. The refactorized functionsS aba b have non-trivial collinear-soft and soft components. The final result is the factorization formula in eq. (4.26) with index j summed over j = 1, 2, 3, 4. In section 5 the hard coefficients H S from tree-level matching will be given explicitly.

Leading log resummation for the endpoint region
Having stated the factorization formula, and discussed the physical intuition that underlies it, this section tackles the resummation of large logarithms of m W /(M χ (1 − z)), m W /(M χ √ 1 − z), and m W /M χ . In section 5.1, we present the one-loop anomalous dimensions obtained by computing the real and virtual corrections to the factorized functions presented in the previous section. We also check consistency conditions for these anomalous dimensions (namely that they sum to zero), thus verifying our factorization formula at the one-loop level. In section 5.2, we describe a simplified resummation path sufficient for LL order and then solve the RGEs and collect all the resummation factors necessary for obtaining the final resummed cross section. The culmination of this work is eq. (5.30). Explicit calculations are given in appendix A. In section 5.3, we demonstrate that our result recovers both the exclusive and inclusive limits.

One-loop anomalous dimensions and factorization consistency
In the results for the anomalous dimensions presented below, we only keep the double log pieces that are required for resummation at LL accuracy. The hard function H i (M χ , µ) only has a µ anomalous dimension, where C A is the SU(2) W quadratic Casimir invariant for the adjoint representation (explicitly C A = 2), i, j = 1, 2, 3, and the structure of the RGE is diagonal. Here, and throughout this section we will useα W = α W /(4π) to simplify the results. Furthermore, we can set z → 1 to leading power, so that the hard function is independent of the infrared measurement. The same anomalous dimension, but derived at the level of the amplitude, was obtained for exclusive heavy WIMP annihilation [42,43,46]. The photon jet function J γ (m W , µ, ν) consists of only virtual diagrams, and is computed in the broken theory. An example diagram is Here the dashed line indicates the final state cut, which puts the single identified photon on shell. We find that the µ and ν anomalous dimensions are given by Due to its fully inclusive nature, we find that it has no anomalous dimension in µ or ν. Instead, these dependences are entirely captured by the matching coefficient H Jn (M χ , 1 − z, µ), which is described by the same diagrams but at the high scale. The dashed line again represents the final state cut, which at NLO can contain one or two particles. Since the one-loop correction to the jet function is a plus distribution, the RG evolution takes a simpler form in Laplace space. We will use s to denote the Laplace variable conjugate to M χ (1 − z). We find its anomalous dimension to be As discussed in section 4.2.3, the general case is complicated by a proliferation of color structures that mix beyond tree-level. For simplicity, we will consider, by top-down construction, only the functions that appear in our analysis at LL order. The µ RGE for thẽ S functions is a matrix equation which exhibits a non-trivial mixing structure. The ν RGE is given by where the matrix is diagonalγS The interpretation of the scales appearing in the functionS = C S S requires some care since this is a combined object. While both the C S and S functions have a natural scale µ = m W (see the ν anomalous dimension given in eq. (5.11)), the scale µ = 1/s appears in the logarithms of the µ anomalous dimension in eq. (5.9). This can be understood from the consistency of the RG, since the µ running of C S and S must combine to yield the natural scale of H S , namely µ = 1/s. Despite its confusing appearance, this appearance of 1/s provides a non-trivial check on our refactorization.

JHEP03(2018)117
One further important feature of the anomalous dimensions in eq. (5.10) is that at LL order, all rapidity anomalous dimensions vanish for µ = m W . We will exploit this feature in section 5.2 by choosing a resummation path where all rapidity evolution is done at the scale µ = m W , eliminating the need for a non-trivial rapidity evolution.
For the matching coefficients H S,ij of the soft sector we have which involves the anomalous dimensions for the soft and jet functions before refactorization, given by As in the case of the hard function, the RG structure for the soft functions S i is diagonal. Using the anomalous dimensions in eqs.
where k, l = 1, 2, 3, 4. One can check that these relations are satisfied using eqs.

Analytic resummation formula
We now have all the necessary ingredients to provide an analytic expression for the resummed spectrum at LL accuracy. As discussed in section 3.2, the resummation can be simplified by making a judicious choice of path in the (µ, ν) plane. Our choice is illustrated in figure 6. Due to the refactorization of the soft function S into the soft and collinear-soft functions, each of which have a complicated color structure, and whose renormalization will involve color mixing, the renormalization group structure is quite complicated for a generic path. However, this can be avoided by noting that at µ = m W , the rapidity anomalous dimensions of the soft and collinear-soft functions given in eq. (5.10) vanish at LL order. Hence, we take the functions at their natural scale -H with µ = M χ , H Jn with µ = 2 M χ /s, and H S with µ = 1/s -and run them all down to µ = m W . Finally, at µ = m W , we can then trivially run the soft, collinear-soft, and jet functions to the same rapidity. This choice of path provides a significant simplification since we can simply compute the µ anomalous dimensions for the functions H, H Jn and H S . Beyond LL accuracy, this is no longer possible, and the full factorization that we have developed in this paper must be utilized.
There is one additional subtlety regarding the evolution structure that has been glossed over in figure 6, but that requires care to reproduce the correct behavior in the limit z → 1. Recall that in deriving our factorization, which is summarized in figure 4, we have assumed the hierarchy In this small region near the endpoint, our EFT is technically speaking invalidated. Physically, the constraint on the final state becomes so restrictive that the jet is composed of a single boson. Due to the intrinsic IR cutoff set by electroweak symmetry breaking, it is unphysical for these scales to go below the scale m W . Instead, we must introduce a Θ-function in the RG evolution, which ensures that the running only contributes in the region where the scales are above m W . As we will see, with this modification, our EFT will correctly transition to the exclusive endpoint calculation. This choice of scales is implemented in (1 − z) space. Therefore, in Laplace space we take arbitrary scales µ H Jn and µ H S (µ H can be set to its canonical value since it is z independent) transform to cumulative space where we can implement our scale setting as a function of (1 − z cut ), and then differentiate to obtain the resummed spectrum. Note that in the following, we will always use z cut when discussing the cumulative space, as per the definition of eq. (1.2).
The RG equations can now be solved in the usual manner. For the hard functions H and H Jn , we derive the evolution kernels where the first and second arguments of the kernels denote the scales we are running between, starting from the natural scale of the relevant function, and ending at µ ∼ m W . For the hard function H S , we need to solve the system of RG equations in eq. (5.13) in order to run from µ = µ H S down to µ = m W . We find that where These kernels resum all leading double logarithms.
To put together the resummed cross section, we need the tree-level values of the hard function H, see eq. (4.12), Here LP −1 denotes the inverse Laplace transform. The prefactors are determined by treelevel matching to full theory, and we have suppressed the arguments of the evolution kernels. At LL accuracy, the cumulative distribution, can be obtained setting s = 1/(2M χ (1 − z cut )) in the Laplace space expression for the cross section, and inserting a 1/(2M χ ) for the measure. At the level of the cumulative, we can now explicitly set our canonical scales as

JHEP03(2018)117
This implements the physical constraint that the jet and soft scales never go below the scale m W . For a more sophisticated analysis, smooth transition functions could be used instead of Θ-functions. This is often done to transition from resummation to fixed order, where the smooth transition functions are referred to as profiles [135]. Here we content ourselves with this simple choice of scales. This simple choice of profiles also allows us to give a closed form analytic result for the differential spectrum involving the Θ-functions. With this choice of scale, the evolution kernels appearing in the cross section, now also explicitly involve the Θ-functions that cut off their evolution as appropriate. For example, for the jet function evolution kernel, we have which becomes unity for m W ≥ 2 M χ √ 1 − z cut . The soft function evolution kernel is completely analogous.
Combining all the ingredients, we arrive at the final expression for the cumulative cross section at LL accuracy Here the Θ-functions explicitly enforce that none of the functions are RG evolved below the scale m W , as emphasized above, and are a crucial part of the final result. Each of the functions appearing in this expression, as well as their physical significance will be defined shortly.
We can now obtain the differential spectrum by taking the derivative of eq. (5.29) with respect to (1 − z cut ). The differentiation of the cumulative result must be performed carefully due to the presence of the Θ-functions, which when differentiated give rise to δ-functions. However, all the δ-functions explicitly cancel, except for the δ-function for the fully exclusive contribution. Carefully performing the differentiation, we obtain the final result for the differential spectrum:

JHEP03(2018)117
This simple formula provides the resummation of all logarithmically enhanced terms to the spectrum at LL accuracy. As before, to simplify the notation, we have written this expression withα W = α W /(4 π). This result is composed of several pieces with clear physical significance, each of which we now explain. The tree-level cross section appears as an overall multiplicative factor, as does the standard massive Sudakov form factor with logarithm The double logarithmic asymptotics is governed by the cusp anomalous dimension [136], in this case at one-loop, where we recall that C A is the Casimir of the adjoint representation of SU (2). Explicitly, in our normalization, C A = 2. In eq. (5.30) we have written Γ 0 in the exponent to emphasize that it is the cusp that controls the anomalous dimensions, but used the explicit form of eq. (5.33) in the prefactors. The first term in the eq. (5.30) is localized at z = 1. Only the Sommerfeld factor |s 0± | 2 appears since the tree-level process is the annihilation of the charged states χ ± . The second term describes the non-trivial z dependence. Here the combination of Sommerfeld factors s 00 s * 0± + s * 00 s 0± , s 00 s * 0± + s * 00 s 0± , (5.34) appear. The perturbative dynamics are controlled by the two logarithms , (5.35) associated with the jet and soft scales, respectively. For convenience, we have also defined Θ functions associated with the range of the soft and collinear scales In addition to the Sudakov logarithms, the z dependence is controlled by the functions

JHEP03(2018)117
which capture the power divergence in 1−z, and the subscript is standard notation denoting that these contain a single power of the logarithm. The presence of the 1/(1−z) factor gives the expected leading power scaling for the cross section. The power divergence for the soft logarithm is cutoff at z = 1 − m W /(2 M χ ) and for the jet logarithm at z = 1 − m 2 W /(2 M χ ) 2 . These physical cutoffs arise from the value of z at which the soft and jet scales hit the scale m W , where the running must be turned off, as has been discussed above. We note that in the massless theory, the power law divergences would be regulated as plus distributions. Instead, here m W explicitly cuts off the divergence at a finite distance from the endpoint.
There is a physical interpretation for each of the different terms in eq. (5.30). The first term, which is localized at the endpoint, corresponds to the fully exclusive cross section, while the other terms describe deviations from the endpoint associated with either soft or collinear radiation. With this understanding of the correct treatment of the scales as we transition to the fully exclusive endpoint, and how they are implemented in our final factorization formula, in the next section we show that our LL expression in the endpoint region correctly reproduces the LL in both the exclusive and OPE regions. Firstly, however, note that expanding eq. (5.30) to fixed order, setting the Sommerfeld factor to its tree-level result |s 00 | 2 = 1, and dropping Θ-functions, we find This result agrees with the O(α 3 W ) logarithm derived in the fixed order calculation of [50].

Reproducing the exclusive and inclusive cross sections
In this section, we demonstrate that our EFT acts as a mother theory which includes both the exclusive (z cut → 1) and inclusive (z cut → 0) results as limiting cases of our resummed expression eq. (5.30). It is important to note that the expansions performed here differ from previous calculations such that power corrections are not expected to be identical. However, this complication is avoided here due to the simple structures that are present at LL order. The focus of this section will be showing how to take these two limits analytically. Section 6 will provide a numerical study of the theoretical error that results from scale variation.

Inclusive limit
To obtain the inclusive limit of the total cross section, we simply integrate the differential cross section given in eq. (5.30) from z = 0 to the endpoint z = 1. Explicitly,

JHEP03(2018)117
Performing the integral, we have where in the last line we have introduced the notation of [41] f This is precisely the result obtained in [41,44,45], demonstrating that we reproduce the inclusive limit to LL order.

Exclusive limit
Note that the signature of interest for experiments like HESS, where the experimental resolution has a width σ m 2 W /(4 M 2 χ ), includes a contribution from the exclusive line and the endpoint spectrum. It is therefore important that we are also able to reproduce the resummed fully exclusive cross section from our factorization. This can be accomplished by integrating eq. (5.30) from z = 1 − m 2 W /(4 M 2 χ ) to z = 1, which corresponds to a kinematic requirement such that only the exclusive final state is possible since both the jet and soft scales are set by the electroweak boson mass. 11 This demonstrates that for the case where the experimental resolution has a width δ m 2 W /(4 M 2 χ ), we have provided the complete description as relevant experimentally (with the additional caveats discussed in appendix B).
When integrating from z = 1 − m 2 W /(4 M 2 χ ) to the endpoint both L J 1 and L S 1 are zero. Therefore, we can we trivially integrate the δ(1 − z) dependent term to find This agrees with the exclusive calculation at leading log accuracy performed in [42,43]. The fact that we reproduce this result makes it straightforward to convolve the resummed photon spectrum with the experimental resolution -no merging between different results is required. In this sense, our EFT acts as a mother theory that completely describes the photon spectrum for heavy WIMP annihilation at LL order. 11 Note that for z > 1−m 2 W /(4 M 2 χ ), eq. (5.30) is proportional to a delta function for exclusive production, namely δ(1−z). It is important to note that we have power expanded away any mass dependence that would lead to kinematic differences between the γ γ and γ Z final states. We therefore are implicitly assuming that the finite resolution function sufficiently smears these differences such that they are not experimentally relevant.

JHEP03(2018)117 6 Numerical results and scale variation
In this section, we provide a numerical study of our final prediction for the spectrum by evaluating eq. (5.30) for wino DM. This allows us to explore the relative contributions from the line annihilation and the endpoint spectrum for different choices of the DM mass. We will also show the cumulative spectra, as given to LL accuracy in eq. (5.29), which provides intuition for the finite bin effects that are relevant to realistic experiments. Then in section 7 we will provide a mock reanalysis of the HESS line search, and will convolve our predicted spectrum with the Gaussian line shape assumed by HESS.
As shown explicitly in section 5.3, our resummed spectrum analytically reproduces the fully exclusive and the fully inclusive limits, so that we can additionally study the transition between these approximations. This clarifies the disparate conclusions that have been drawn using these different approaches. In particular, the exclusive calculations of [42,43,46] claimed a reduction factor of ∼ 2.2 when compared with the tree-level cross section for a 3 TeV wino. For contrast, the inclusive calculation of [41,44,45] found a reduction of only a few percent. Physically, this results from the fact that an increasingly exclusive constraint on the final state implies there will be less cancellation between the virtual and real corrections (for discussions in the context of electroweak logarithms, see e.g. [137,138]). The proper interpretation of the experimental limits depends on how rapidly the transition between the exclusive and inclusive cross sections occurs. Our EFT analysis provides a complete and decisive resolution of this issue. Interestingly, we find that the experimental values of current interest to the HESS line search, z cut ∼ 0.8-0.9, lies right in a transition region between the two limiting cases. This emphasizes the need to properly treat the impact of finite resolution, as we will do in the next section. However, before moving to our mock reanalysis of the HESS search, we will provide some numerical results along with an estimate of the impact of scale uncertainty.
In figure 7 we show the differential spectrum z 2 d σv /dz for several values of the DM mass. The delta function contribution from the exclusive process is not included. We see that the endpoint tracks the mass of the DM as expected. Furthermore, the contribution from the resummed continuum grows as the DM mass is increased. However, this effect is mitigated by the strong mass dependence of the overall cross section, both due to Sommerfeld enhancement and the overall 1/M 2 χ scaling, which explains why the 3 TeV result lies above both the 1 TeV and 10 TeV results. The kink in the 1 TeV distribution is a result of the Θ-functions appearing in the choice of scales, as discussed in section 5.2 (in reality, there are kinks in all the distributions, but they are only visible by eye for the 1 TeV distribution). This kink is ultimately unphysical and could be removed by a smooth choice of scales, but is well within our uncertainty bands.
The uncertainty bands in figure 7 are the result of varying the renormalization scales corresponding to the natural scales of the functions appearing in our factorization. Due to our choice of renormalization path, we simply vary the µ scale of the different functions by a factor of two about their natural scales.
An alternative numerical representation of our results is provided in figure 8, where we plot the cumulative cross section as a function of the z cut , for several values of the DM mass. Here we do include the delta function contribution that yields the exclusive annihilation process, which accounts for the finite value when z cut = 1. The uncertainty bands are computed using the same prescription for the scale variation performed for figure 7.

JHEP03(2018)117
The two endpoints, namely z cut = 1 and z cut = 0, correspond to the fully exclusive and fully inclusive limits, respectively. Interestingly, for the experimentally relevant range z cut ∼ 0.8-0.9, the cumulative cross section takes an intermediate value approximately midway between the two extremes. This implies that for these values of z cut , logarithms of the resolution are playing an important role, in keeping with the conclusions of the fixed order calculation in the inclusive limit [45]. Theoretically robust results require the all-orders resummation of logarithms from finite bin effects, as has been done here for the first time.

Impact on indirect detection constraints
The resummed photon spectra derived above have clear implications for heavy DM line searches. In particular, thermal wino annihilations would produce TeV scale photons. 12 When these photons strike the Earth's atmosphere, they initiate a detectable shower of particles that persists to the surface. Exactly reconstructing the energy of the incident photon from the resultant shower is impossible, and as such any real instrument will need to account for finite energy resolution effects associated with the spread of possible reconstructed energies given a single true energy.
As discussed at the outset, the strongest constraints on the wino are due to HESS observations of the Galactic Center [26,28]; updated limits are expected shortly involving the full HESS I dataset [139,140]. Line searches are typically designed to be modelindependent, and thus assume that only the line emission is relevant (although some specific non-line hard spectra have also been tested [28,141]). As we demonstrated in figure 8 above, photons away from the endpoint contribute to a finite bin at a non-trivial level. This is especially true for HESS, where the effective z cut ∼ 0.8-0.9 depending on the incident energy. Furthermore, the line analysis of HESS is not a bin-based counting experiment but requires subtraction of an unknown background, which is modeled by a smooth function. The presence of signal photons at even lower energies may bias the data-driven background model if this signal spectrum is not correctly modeled, further modifying the limit.
The goal of this section is to estimate how much including the correct shape and normalization of the resummed spectrum would be expected to change the HESS limit, relative to the case of a pure line.
It is important to emphasize that the results presented in this section are approximate, and should not be taken as updated limits on the wino. At issue is that the full dataset HESS used to construct their limits in ref. [28] is not public. What we will show are results from a simplified mock version of that analysis, using a Gaussian likelihood rather than the full likelihood, which has been validated to yield comparable limits when assuming exclusive line emission. We can then explore how the various conclusions are modified when we include the endpoint emission spectrum. The conclusion is that the additional emission should strengthen limits on the wino by a O(1) factor. This provides motivation for future experimental analyses to include these contributions when determining limits.

JHEP03(2018)117
This section contains three parts. First, we review how to map from DM model parameters, including the relevant astrophysical inputs, to a prediction for the number of photons that HESS would observe. Then we apply this formalism to demonstrate the range of parameters that HESS can constrain. Finally, we outline our mock analysis procedure and present approximate results showing the impact of our resummed spectra on current constraints.

Predicting the indirect detection flux
In order to determine the sensitivity to wino DM, we need a prediction for the number of photons that should arrive at an experiment as a function of the DM parameters. This can be derived using the canonical indirect detection formula, which specifies the differential energy flux arriving at the detector, where Ω ROI ≡ ROI dΩ. The particle physics contribution σv /(8 π M 2 χ ) dN γ /dE depends on the velocity averaged total annihilation cross section σv , which is summed over all final states involving a photon, and the average photon spectrum per annihilation, dN γ /dE, which can be written as 13 dN where the f index refers to the different final states with associated branching fractions Br f and photon spectra dN f γ /dE. Since the spectrum here is the result of resumming multiple electroweak final states (not including the photons that result from decay of unstable W ± and Z bosons, see appendix B for a discussion), we will only refer to the total averaged quantity dN γ /dE for the remainder of this section.
The remaining ingredient is the so-called J-factor, which is an astrophysical input. It is determined by the distribution of the DM along the line of sight in the region of interest (ROI) under consideration. It additionally accounts for the fact that two particles must find each other for annihilation to occur; the J-factor depends on the number density squared as where ρ DM is the Milky Way DM mass distribution, s is the distance from Earth along the line of sight, and Ω gives the coordinates on the sky within the ROI. Note that as written, the J-factor has units of TeV 2 · cm −5 , and in particular there are no units of sr due to the denominator in eq. (7.3). We caution, however, that a number of other conventions are in use. 14 For a fixed ROI, J is then in principle determined by the Milky Way DM profile. Unfortunately, the shape of ρ DM is very uncertain near the Galactic Center, see e.g. [143],

JHEP03(2018)117
and in particular within the ROI of the HESS search of ref. [28]. For the case of the wino, once the mass is fixed the cross section is fully specified. Therefore, one can translate limits on wino annihilations into a constraint on J, as done in figure 10 below.
It is also of interest to fix a prototypical value for J and then set a limit on the annihilation cross section, since this is how these constraints are typically presented. For this purpose we adopt the Einasto profile, the default profile assumed in the HESS analyses, which is given by where r is the distance from the center of the halo, and following [144], by default we take α = 0.17, r s = 20 kpc, and then normalise the profile so that we reproduce the local DM density of 0.39 GeV cm −3 at our location which is 8.5 kpc from the Galactic Center. Another frequently invoked DM distribution is the Navarro-Frenk-White (NFW) profile [145], which takes the form where again we take r s = 20 kpc. We will also make use of the NFW profile (including the possibility of a non-trivial core) when interpreting our results below. Finally, putting this all together results in the differential energy flux arriving at the detector, 1 Ω ROI dΦγ dE , which has units of photons · cm −2 · s −1 · TeV −1 · sr −1 . This quantity can be converted into a predicted number of photons (per unit area per unit time) arriving at the experiment from DM annihilation by first multiplying by the solid angle of the considered ROI, Ω ROI , and then integrating over the energy range determined by the experimental search. This photon flux Φ γ has units of photons · cm −2 · s −1 . Converting this to the actual number of photons depends on the experimental effective area and time over which the ROI is observed; a larger detector and longer observations will result in more observed photons. For HESS, the effective area is ∼ 10 9 cm 2 at 1 TeV and current searches make use of 112 hours of observations of the Galactic Center, yielding sensitivity to fluxes ∼ 10 −14 cm −2 s −1 . We can then constrain the DM model using this prediction for the number of photons as an input to a likelihood analysis.

From predictions to constraints
Before we give the details of and results from our mock analysis, it is useful to discuss how we are mapping from the theory prediction to the experimental constraints. The subtlety arises because the original search was performed under the assumption that the annihilation signature is a line; in this case, by definition all photons have the same energy. The spectrum of a typical WIMP can be decomposed into two contributions The line is due to exclusive annihilations to γ γ and γ Z. Since our interest here is in heavy WIMPs, we will neglect the fact that the finite Z mass causes for the photons that result from the γ Z process, and will combine these line contributions using σv line ≡ σv γγ + 1 2 σv γZ , (7.7) with E γ = M χ for all line photons. The continuum receives many contributions. In the DM literature, this is usually separated into photons from "internal bremsstrahlung" [47][48][49][50][51], as well as final and initial state radiation, on one hand, and those photons that result from the cascade decay chain of unstable particles on the other hand. The decay processes can yield many final state photons with a broad energy spectrum. Our endpoint calculation for winos resums the nondecay perturbative processes, and as such it does not include the additional continuum photons that result from the decay of the W ± and Z. However, this contribution is demonstrated to have little impact on the limits for heavy winos in appendix B. This conclusion is intuitive since the photons from the W ± /Z cascade decays are much lower energy than the exclusive and endpoint contributions. Therefore, we model the continuum as only being due to the endpoint contributions, which we denote with E(E), and using eq. (5.30) the LL result is given explicitly by where as usual, z = E/M χ . The resulting spectrum per annihilation is such that σv line / σv is the branching fraction to line photons. Note that our calculation predicts not only the shape of the endpoint contribution, but also the relative normalization of this with respect to the line spectrum. Putting these details together, we arrive at the theory prediction which is idealized in the sense that it neglects experimental effects. As such we are still missing one ingredient, which is the fact that we need to convolve this with the experimental energy resolution. We can describe the energy resolution via a convolution function Σ(E − E ), where E is the true photon energy and E is the reconstructed value, and the spectrum an experiment would measure is

JHEP03(2018)117
The HESS collaboration has published a model for Σ(E −E ) which we use here, a Gaussian that is peaked near the true energy with a width that varies from 17% at 0.5 TeV and 11% at 10 TeV. We interpolate in between these values using the log of the energy, and find a width ∼15% at 3 TeV. HESS can constrain the overall normalization of eq. (7.11); in terms of the theory prediction, this can be interpreted as a constraint on the quantity However, it is critical to specify the assumed energy spectrum E(E) (in addition to a line contribution) when deriving a HESS constraint on the cross section. For the following comparisons, we will use the LL endpoint spectrum computed in this work, so that where C HESS is the coefficient that is constrained using the HESS data, we take E LL E from eq. (7.8), and Σ(E − E) is as discussed above. In the next section, we will interpret the HESS data as a constraint on C HESS using a mock analysis, and will then convert this into an approximate constraint on winos using eq. (7.12). We will either use eq. (5.42) to predict σv line for a given mass in order to set a constraint on J, or we will assume the Einasto profile which gives us J and then constrain the cross section σv line . We will also provide a constraint on the core size, using the NFW profile modified to include a core. Note that we can test the effects of ignoring the non-line endpoint contributions by simply setting E(E) = 0; up to the approximations in our analysis required by not having the full likelihood available, this should reproduce the limits stated in ref. [28]. This allows us to directly compare constraints on the line only and the line plus endpoint spectrum, thereby highlighting the impact of our main result eq. (5.30). The next section outlines the details of our mock analysis and provides approximate constraints on either the cross section or the J-factor.

Approximate constraints
Using the procedure described in the previous section, one can in principle interpret the HESS data as a constraint on wino DM annihilations. As the data collected by the instrument is not public, we are not able to provide a full and precise update of the constraints on winos. Instead, we will perform a simplified mock version of the HESS analysis in order to estimate the impact of the corrections calculated here on the resulting limits. Our mock analysis can roughly recover the published line limits in the case where we take E(E) = 0 above. We will then extend the analysis to include the endpoint contributions, demonstrating that they strengthen the limits by an O(1) factor.
Our mock analysis is based on a simplified version of the analysis performed in ref. [28]. Figure 1 of that work provides the measured flux and the associated uncertainty as a function of energy in their ROI near the Galactic Center. We digitized this dataset and used it as the input to a Gaussian likelihood analysis. We note that since HESS is an JHEP03(2018)117 instrument that counts the number of incident photons, the Poisson likelihood should in principle be used. However, the number of counts cannot be exactly reconstructed from the publicly released flux data. The non-Gaussian nature of this dataset is made manifest by the asymmetric error bars that are particularly clear at higher energies. We approximately included the asymmetry in the likelihood by using the upper error bars to determine the likelihood contribution from bins where our model prediction exceeded the data, and the lower error bars for bins where our model prediction fell below the data. We found that this approach gave better agreement with the published HESS results than symmetrizing the error bars.
The dataset d i is defined using energy bins with associated index i, where the digitized HESS flux gives a central value µ i and error σ i , chosen (between the upper and lower error bars) in the manner described above. The prediction m i (θ) is a function of the model parameters θ. The DM-signal contribution to the model is computed using eq. (7.11). We will treat this theory flux as being a function of the DM mass, M χ , the line photon cross section, σv line , and the J-factor. As emphasized above, given M χ we can either calculate σv line and then constrain J, or assume a value of J and turn this into a constraint on σv line .
Even in the most optimistic DM scenario, the events collected by HESS will not be solely due to DM annihilation. Firstly, there is a substantial flux of cosmic rays colliding with the atmosphere, which can mimic gamma-ray signals. Secondly, there will be genuine gamma-rays due to high-energy astrophysical processes, such as protons in the inner galaxy colliding with gas and producing energetic neutral pions which decay to gamma-rays. The expected flux from cosmic-rays and astrophysical sources of gamma-rays is not well understood in the HESS energy range, and as such ref. [28] parametrized the background contribution using the following seven parameter model: (7.14) The background is then specified by the seven parameters θ bkg = {a 0 , a 1 , a 2 , a 3 , β, µ x , σ x }.
Combining the signal and background, we arrive at our full model prediction of so that the model is specified by three signal and seven background parameters. From here, given the HESS dataset described above, d = {d i } = {µ i , σ i }, we can write down our assumed Gaussian likelihood function as In order to restrict our likelihood to be a function of only the signal parameters, we eliminate the nuisance parameters using the profile likelihood method, L d|θ sig = L d|θ sig ,θ bkg , (7.17) where the hat indicates evaluating the function at the values of θ bkg that maximize the likelihood (see [146] for a review).
With this reduced likelihood, we can then define a test statistic for upper limits as a function of M χ on either σv line or J. To begin with, we can fix J and set a limit on σv line . To determine the fixed value of J, we use eq. (7.3) assuming an Einasto profile as given in eq. (7.4). The ROI for this dataset was a 1 • circle around the Galactic Center, with the Galactic plane masked for latitudes less than 0.3 • , which yields J 7.39 × 10 18 TeV 2 cm −5 . (7.18) Fixing this value, we define the test statistic as where again a hat denotes the value that maximizes the likelihood. Using this test statistic, the 95% limit on σv line is then determined by solving for q σv line (M χ ) = −2.71, and is shown in figure 9. In this figure we have also shown the prediction for the wino cross sectionif these were exact limits and if the DM distribution followed an Einasto profile in the inner galaxy, then the wino would be excluded when this prediction is above the mock limit curve.
This figure also contains the published HESS limits, taken from figure 4 of [28]. The extent to which our line-only limits disagree with the published values highlights that our mock analysis is not exact and thus should not be taken as the true limit on wino DM. Nevertheless the figure does make it clear that the addition of the endpoint contributions can lead to a non-trivial enhancement on the sensitivity. For this reason, the effects calculated in this work represent an important contribution that should be included in future searches for heavy DM annihilation.
Alternatively, for limits on J, we fix σv line to the exclusive wino prediction appropriate for that mass using eq. (5.30), and in a similar notation to [147], define our test statistic as As for the cross section, this test statistic allows us to establish the 95% limit at a given mass through the relation q J (M χ ) = −2.71, and the result is shown in figure 10. In this case we have repeated the analysis with and without the endpoint contributions calculated in this work, with the impact of our calculation being as much as a factor of 3 improvement in the limit, and a factor of ∼1.5 at the thermal mass.  Figure 11. The NFW core size required to save the wino as derived using our mock analysis. This figure follows from the J-factor constraints given in figure 10. At a given mass, the constraint on J can be converted into a core size limit by calculating the corresponding cored NFW J-factor in the HESS analysis ROI. For a thermal wino at 3 TeV the estimated constraints improve from 0.70 to 0.99 kpc when including the endpoint contributions. This core size is beginning to be probed in both numerical and astrophysical settings. We again emphasize that these constraints should be only taken as an estimate.
The results above demonstrate that updating the wino limits to include the endpoint contribution can easily lead to O(1) improvements in the limits on σv line or the Galactic Center J-factor. Finally, we emphasize that the search for the wino is reaching a level of sensitivity such that O(1) factors are important. One way to see this, is by converting the limits into a statement on how large a core in the Milky Way DM density profile is required to save the wino from the HESS constraints.
For concreteness, we use a cored version of the NFW profile, following [24]. For r > r core , we take the NFW profile as defined in eq. (7.5); when r ≤ r core , we set the profile to a constant value ρ NFW (r core ), such that the density profile is flat within the core radius. Restricting ourselves to cores smaller than 8.5 kpc, the presence of a core reduces the associated J-factor of the halo. In this way we can directly convert J-factor limits into a corresponding constraint on r core , which we show in figure 11. From this, we can see that for a thermal wino at exactly 3 TeV, the estimated core constraint increases from 0.70 kpc to 0.99 kpc when including the additional photons from the endpoint spectrum.
The values constrained in figure 11 turn out to be at the edge of the core sizes that are beginning to be probed by a combination of numerical simulations and data. On the numerical side, it was shown that recent simulations of Milky Way-like halos in simulations including the effects of baryons, can potentially contain cores up to O(1) kpc [148]. The total DM mass in the Galactic Bulge region can be estimated from observations of stars JHEP03(2018)117 in the Bulge [149], and disfavors a canonical NFW profile with a core size larger than ∼2 kpc [150]. The core sizes needed for the thermal relic wino to survive indirect detection bounds are thus beginning to be constrained by stellar observations; accounting for the detailed endpoint spectrum is an important component when drawing this conclusion.

Conclusions
In this paper we have developed a comprehensive effective field theory framework to compute the photon spectrum for annihilating (or decaying) DM. We provided a new factorization formula, which allows for a resummation of all large logarithmic contributions, properly treating the effects due to electroweak symmetry breaking, the experimental resolution on the γ + X final state, and the Sommerfeld enhancement. We have computed the relevant one-loop anomalous dimensions, showing the consistency of the factorization formula at this order. We have shown that the contribution from the spectrum has a numerically important effect for experimental searches of interest, e.g. gamma-ray line searches from the HESS telescope. Our final result is a compact analytic expression for the differential spectrum at LL accuracy, which can easily be convolved with experimental resolution functions to provide realistic predictions.
Our EFT can be interpreted as a mother theory that includes as particular limits the fully exclusive and fully inclusive cases. The framework presented here correctly describes the transition between these two limits, allowing us to understand how Sudakov double logarithms impact the spectrum as a function of the experimental resolution. It also allows us to assess the extent to which a fully exclusive or fully inclusive approximation, as had been previously considered in the literature, is appropriate. Interestingly, we find that for the range of resolution parameters applicable for current and near future experiments, the result is intermediate between the fully exclusive and fully inclusive predictions. This resolves the differing conclusions obtained in the literature, and provides a unifying picture of the importance of Sudakov resummation for indirect detection searches. We have estimated the impact on the interpretation of current searches by providing a mock reanalysis of the HESS data, and we find that we are probing core sizes in a region where precise calculations of the particle physics components are relevant. Now that this paper has established an EFT framework for describing the photon spectrum resulting from DM annihilation, one can extend this work in a number of future directions. It would be of formal interest to understand the structure of the factorization and resummation at higher logarithmic order. Although the electroweak couplings are small, significantly improved uncertainties have been observed at NLL [42,43,46], implying that NLL is likely the highest order that is relevant. Additionally, the explicit NLO calculations provided in [42,46] demonstrate that higher order terms that are not logarithmically enhanced are numerically unimportant, justifying our choice to neglect them.
There are also additional phenomenological applications. One could extend these results to other heavy DM models, e.g. a thermal Higgsino, a mixed wino-higgsino, or minimal DM candidates. In many of these cases, the constraints can be different [24,151,152], implying that a dedicated analysis is warranted. From the point of view of extending the JHEP03(2018)117 work presented here, the underlying EFT is unchanged, but one must modify the Sommerfeld calculation and the explicit values for the hard matching coefficients and anomalous dimensions.
A simple heavy DM candidate provides a viable and phenomenologically relevant explanation for the observed relic abundance that could show up in current or future indirect detection searches. This work casts the prediction for the photon spectrum that can result from this class of models in a theoretically robust setting, where perturbation theory can be maintained by performing resummation of all large double logarithms. If a signal of heavy DM annihilation appears, this work will be critical to interpreting it.
Photon jet function. The photon jet function, which has a single photon as its final state, is defined in eq. (4.13) as Evaluating this function at one-loop yields where µ and ν are the virtuality and rapidity renormalization scales respectively. Here n f denotes the number of fermion flavors. We take n f = 5 in our numerical results. The µ and ν anomalous dimensions can immediately be extracted from this result, and we find, Recoiling jet function. When computing the recoiling jet function, all IR divergences are explicitly regulated by the measurement of the final state mass. This is unlike the photon jet function, where the scale m W acts as a regulator. To compute the anomalous dimensions, it is therefore sufficient to expand away the scale m W from the beginning, simplifying the calculation. The inclusive recoiling jet function is then defined as We can rewrite this jet function as

JHEP03(2018)117
with p = (2 M χ , k + , 0) in order to enforce the δ-function measurement constraints. Written in this form, the function is completely inclusive. Therefore, we can use the optical theorem to write this as the imaginary part of the forward scattering amplitude This jet function has been evaluated in the literature [153][154][155]; the one-loop result is where the subscript plus denotes a plus distribution, see e.g. [156] for an extensive review of its properties. The kinematics for heavy DM annihilation imply that p 2 = 2 M χ k + , so that To expose the simple renormalization group structure, we transform to Laplace space, where the Laplace conjugate variable of k + is taken to be s. Keeping only the leading log term, we find where γ E is the Euler-Mascheroni constant. Finally, we extract the µ anomalous dimension This function has no rapidity anomalous dimension as it is a SCET I type function.
Ultrasoft function. There are four operators that contribute to the ultrasoft function in the EFT: S 12 , S 21 , S 11 , S 22 , see eq. (4.11) above. Using the expressions below, we can then extract the LL µ and ν anomalous dimensions. We will find that each operator yields the same result, This calculation will also expose additional IR divergent contributions, which is the sign that refactorization is necessary. The one-loop results will be expressed in terms of several integrals, denoted in bold and labeled with V and R for virtual and real respectively, which are defined and evaluated below. These integrals are evaluated using dimensional regularization as an IR regulator,

JHEP03(2018)117
and with the rapidity regulator as defined in section 3.2. The integrals that we will require in our calculation are defined as follows. The nn integrals are the nv integrals are and thenv integrals are Next, we consider each of the four ultrasoft functions in turn. First we provide the operator definition, followed by the tree-level and one-loop evaluation in order to compute the anomalous dimensions for the different color structures of the ultrasoft function. Since we are doing this in the EFT before refactorization, we will refer to these as ultrasoft functions and the corresponding states as |X U S . These ultrasoft functions will ultimately be refactorized.
• S 11 is defined as

JHEP03(2018)117
Evaluating at tree-level in Laplace space yields S aba b 11 tree = δ ab δ a b , (A. 23) and at one-loop yields where the second line is expressed in Laplace space. Extracting the LL anomalous dimensions from these results yields eq. (A.15).
• S 12 and S 21 are defined as Evaluating at tree-level in Laplace space yields and at one-loop yields where the second line is expressed in Laplace space. Extracting the LL anomalous dimensions from these results yields eq. (A.15).
This result manifests the same UV virtuality and rapidity divergences as in the case of the S 11 operator which is why it yields the same anomalous dimension as S 11 . However, we see an additional IR divergence appears in the form of log 2 m W s . This results from the non-singlet nature of this operator. In order to factorize this new double log, we need to match this ultrasoft operator onto an EFT below the scale s. This allows us to separate the scales s and m W , yielding our final fully factorized result.
• S 22 is defined as

JHEP03(2018)117
Evaluating at tree-level in Laplace space yields S aba b 22 tree = δ a3 δ a 3 δ bb , (A. 29) and at one-loop yields where the second line is expressed in Laplace space. Extracting the LL anomalous dimensions from these results yields eq. (A.15). Note that although the result appears not to be symmetric in the color structure, the wavefunction F a b ab defined in eq. (5.24) is symmetric under the interchange a, a ↔ b, b .

JHEP03(2018)117
which is given in eq. (A.14), and d d log µ Jn m W , µ, ν = 0 . (A. 34) To the order that we work, we need just the tree level value for Jn, which is Anomalous dimensions for the refactorized ultrasoft function. Unlike for the jet functions, the refactorization of the ultrasoft function is significantly more involved. As given in the text, the general form of the refactorization is .
The goal of this section will be to describe this refactorization in more details, and derive the required anomalous dimensions. Before considering the structure of the anomalous dimensions, we must first derive the color structures of the collinear-soft and soft functions, which were only stated without derivation in the main body of the text. The structure of the Wilson lines in the soft and collinear-soft functions can be derived by performing the BPS field redefinition iteratively. We therefore return to the two amplitude level operators (see eq. (4.2) above) Next we iteratively perform the BPS field redefinition for both the collinear-soft modes and refactorized soft modes, We can now derive the soft and collinear-soft functions in the standard way, by squaring the amplitude level operators and settingD =D ,C =C = 3. We find To simplify the notation and focus on the color structures, we have suppressed the measurement function.

JHEP03(2018)117
One additional complication that arises in the refactorization of the ultrasoft function, is that there are non-trivial zero-bins [106] that must be correctly incorporated. We therefore briefly discuss the structure of these zero-bins, and their dependence on our choice of regulator, showing through two examples how the factorization correctly reproduces the structure of integrands once the zero bin is included. We consider one example of a virtual integral and one example of a real integral, arising from the S 11 integrand.
• Consider the virtual integral (see eq. (A.17) for the evaluation of the unexpanded integral) Let us now consider the collinear-soft limit ( + − ) of this integral. It would appear that according to the power counting the only effect is to drop − from the rapidity regulator term | + − − | η/2 . Since the rest of the integrand is unchanged, this would lead to an unregulated divergence as − → ∞. We would then be forced to introduce a new regulator to counter this divergence. While there are several ways to do this (a ∆-regulator [157], for instance), the simplest way is just to keep the original form of the rapidity regulator. The choice of the regulator we use will affect the zero-bin subtraction that will be needed. so that we recover the full US virtual contribution.
• Now, consider the real emission integral (see eq. (A.16) for the evaluation of the unexpanded integral) (A.42) The soft limit is The CS limit is identical to the full integral. Applying the zero-bin subtraction to this (which turns out to be the same as the soft integral), we are left with Thus, once again we recover the full US real contribution adding together the CS and soft limits. The zero-bin subtraction scheme then is to simply subtract the soft limit of the CS integrals from the CS sector.

JHEP03(2018)117
The analysis of these integrals provides a non-trivial check that our factorization is indeed correct, and that the infrared is completely reproduced by our factorized description.
Having understood the operator basis and the structure of the zero-bin subtractions, we can now compute the anomalous dimensions of the functions arising after the factorization of the ultrasoft function. Here we can considerably simplify the calculation by using the choice of resummation path described in section 5.2 and shown in figure 6. In particular, for this path it is not necessary to separately run the collinear-soft and soft functions. We can therefore simplify our refactorization to and only compute the anomalous dimensions for the functions H S,ijkl and S kl aba b . This drastically simplifies the calculation, since the structure of the color mixing for the collinearsoft and soft operators is quite involved. In the remainder of this appendix we give the explicit results for the anomalous dimensions for H S,ijkl and S kl aba b for all relevant color channels appearing in our factorization.
For ease of notation, as in the body of the text, we will define our ultrasoft operators as, see eq. (4.12), (A.46) In this notation, the refactorization of the ultrasoft function is given by, see eq. The tree-level, and one-loop results, along with the µ and ν anomalous dimensions for the H andS functions appearing in the factorization are as follows: •S 1 is defined as × X c S δ q + − P + X 3f n V dg n (0) 0 δ f g δ a b δ f g δ ab , (A. 48) where the soft sector Wilson lines have contracted to the identity. By inspection, the anomalous dimension for this operator is identical to S aba b 1 , implying that H S,11 = 1 is the only non-zero matching coefficient.
•S 3 is defined as The second line here is essentially the IR piece of the term log 2 m W s . Extracting the LL anomalous dimensions yields d d log µS which shows a mixing betweenS 3 andS 1 , along with d d log νS (A.53) We can now read off the matching coefficients + δ bb δ aa − 3 δ a 3 δ a3 α W π −2 log µ m W log µ s + log 2 µ m W .

JHEP03(2018)117
B Impact of continuum photons from cascade decays In the main body of this work we presented a calculation of the internal bremsstrahlung (+ initial/final state radiation), or endpoint, contribution to the wino annihilation spectrum. As we mentioned there, another source of photons arises from the final state decay products of the unstable particles that are produced by DM annihilations, such as the W ± and Z bosons. In this appendix we estimate the contribution from these additional final states, and show that they have a small impact on the HESS constraints for the thermal wino. However, they could be interesting for instruments searching for lower energy photons such as Fermi.
In order to estimate these contributions, we have added to the line and endpoint spectra the spectrum coming from the decay of W ± and Z bosons. The spectrum of photons that arises from their decay is determined using PPPC4DMID [158] with electroweak corrections turned off, 15 whereas the branching fraction is evaluated differently for the two cases. For annihilation to W + W − , the branching fraction is given by the Sommerfeldenhanced tree-level cross section for this final state [3,4]. 16 Radiative corrections to this cross section, which have been shown to be small [40], are not included. To determine the Z production cross section, we use the leading log cross section, which is given by eq. (5.42) reweighted by c 2 W /s 2 W . In figure 12, we show the impact on the photon spectrum from DM, after convolving it with the HESS energy resolution, when this continuum contribution is added, for two DM masses. Generically, as we approach E γ ∼ M χ , this continuum emission is a sub-dominant effect. However, at lower energies it can have substantial impact (note this spectrum is multiplied by E 2.7 which downweights the flux at lower energies). Nevertheless, such a contribution over many energy bins is hard to distinguish from the 7 parameter background model used by HESS. These background parameters are profiled over, so that we would not expect this additional emission to make a sizable impact. Indeed, in figure 13 we demonstrate this point, by repeating the analysis from section 7.3 with the inclusion of the additional continuum photons. We note the effect of including the continuum becomes more important at higher masses, but is almost always subdominant to the impact of adding in the endpoint emission. Further, the broad nature of the continuum emission can lead to a non-trivial interplay with the background model in fits to the data, and in fact lead to weaker limits for some masses. For example, near M χ = 9 TeV in figure 13, the Sommerfeld resonance leads to an enhancement in the continuum emission. This additional emission at lower energies drives down the best fit background model, resulting in a reduced background prediction near the dark matter mass where the line and endpoint contributions dominate, and accordingly a weaker limit. 15 The electroweak corrections in PPPC4DMID include a partial accounting of the endpoint corrections that we determined in the main body, which they include following [159], and so we remove them to avoid double counting. We thank Marco Cirelli for confirming this point. This choice means we are missing the electroweak corrections from the remainder of the continuum spectrum, however we have confirmed these effects are small by directly comparing the spectra to the predictions of Pythia 8.219 [160][161][162], which includes electroweak showering [163]. 16 Note that there is a factor of 2 missing in the off-diagonal terms of Γ W + W − in eq. (28) in ref. [3], which is corrected in ref. [4].  Figure 13. As in figures 9 and 10, but showing the impact of adding the continuum contribution from W and Z decays in addition to the endpoint on the constraints. In general these contributions have a much smaller impact than that already resulting from adding in the endpoint spectrum. The exception is near the Sommerfeld resonances, where the associated enhanced continuum emission is imprinted on the limits. We caution once more that these are only estimated limits.

JHEP03(2018)117
Finally we note in passing that the large contribution from the continuum may be relevant to lower energy instruments such as Fermi -LAT. The advantage of such an approach is that we can look at a number of different potential astrophysical sources of DM flux, each associated with partially uncorrelated systematics on their J-factors. In this way we can extend the search beyond the Galactic Center and its large uncertainties to look at potentially cleaner environments such as the Milky Way Dwarfs [164,165] or even galaxy clusters [142,166]. However, note that the effective area of Fermi -LAT drops sharply at TeV energies. This implies that if the DM mass is multi-TeV, the HESS constraints are generally stronger than those from Fermi. Even accounting for the astrophysical uncertainties, the HESS dataset continues to be the best probe of the gamma-rays from thermal wino DM.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.