Precision Photon Spectra for Wino Annihilation

We provide precise predictions for the hard photon spectrum resulting from neutral SU$(2)_W$ triplet (wino) dark matter annihilation. Our calculation is performed utilizing an effective field theory expansion around the endpoint region where the photon energy is near the wino mass. This has direct relevance to line searches at indirect detection experiments. We compute the spectrum at next-to-leading logarithmic (NLL) accuracy within the framework established by a factorization formula derived previously by our collaboration. This allows simultaneous resummation of large Sudakov logarithms (arising from a restricted final state) and Sommerfeld effects. Resummation at NLL accuracy shows good convergence of the perturbative series due to the smallness of the electroweak coupling constant - scale variation yields uncertainties on our NLL prediction at the level of $5\%$. We highlight a number of interesting field theory effects that appear at NLL associated with the presence of electroweak symmetry breaking, which should have more general applicability. We also study the importance of using the full spectrum as compared with a single endpoint bin approximation when computing experimental limits. Our calculation provides a state of the art prediction for the hard photon spectrum that can be easily generalized to other DM candidates, allowing for the robust interpretation of data collected by current and future indirect detection experiments.


Introduction
Indirect detection is critical to the hunt for multi-TeV WIMP dark matter (DM). New data are continually being collected by current experiments, e.g. H.E.S.S. [1][2][3], HAWC [4][5][6], VER-ITAS [7][8][9], and MAGIC [10,11], and a number of dedicated line searches for photons have been performed [3,12]. Future experiments such as CTA [13,14] will provide even greater sensitivity. Deriving the experimental ramifications these data will have on the parameter space of DM models requires precise predictions for the hard photon spectrum. Due to finite resolution effects inherent to the relevant experiments, a reliable prediction for not only the rate but also the shape of the spectrum is required to derive robust comparisons between theory and experiment [15].  TeV [32], as a function of resolution parameter z cut , showing our results at both LL and NLL. The NLL calculation significantly reduces uncertainties as compared to LL. A region appropriate for the H.E.S.S. experimental resolution, which is ∼ 10% at these energies [33], is shown in the grey band. This band is representative of the range of values that will contribute when our spectrum is convolved with the H.E.S.S. energy resolution function.
In [15] we extended these EFT approaches to allow for the calculation of the hard photon spectrum in the endpoint region, where the photon energy E γ is near the DM mass M χ , as is relevant for line searches. Our framework additionally allows for the resummation of resolution effects α n W log m 1 − z with m ≤ 2 n − 1, where z = E γ /M χ . Such logarithms are directly related to the experimental energy resolution, since z quantifies the distance from the exclusive case, given by a line at z = 1. A finite experimental resolution smears photons with a small 1 − z into the expected exclusive event rate, and our calculation is able to realistically incorporate such effects.
In this paper we extend this result to next-to-leading logarithmic (NLL) accuracy. To understand the importance of including the NLL corrections, we note that for the situation of interest here the logarithms L become large enough so that L 2 ∼ 1/α W . A leading logarithmic (LL) calculation then captures all terms scaling as 1, and so should provide a good description of the shape of the distribution. However, a LL calculation does not probe higher order radiative corrections, and therefore typically has large uncertainties. On the other hand, an NLL calculation captures the first radiative corrections scaling like α W L, and therefore typically provides a large reduction of the theoretical uncertainties. This reduction of uncertainties is clearly illustrated in Fig. 1, which shows a comparison of our earlier LL calculation with the NLL result achieved here. With NLL accuracy, the theory uncertainties become a subdominant contribution to the total uncertainty relevant experimentally.
While our calculational framework is generally applicable to any heavy WIMP candidate, we will often specify to the case where the DM candidate is a wino, the neutral component of a triplet of SU(2) W with zero hypercharge. In addition to allowing us to illustrate our formalism, the wino is well motivated phenomenologically (see e.g. [34][35][36][37][38][39][40][41][42][43]), making our results of interest to current experiments. In a companion paper [44] we have used these results in a realistic H.E.S.S. forecast analysis to study the impact of having a complete description of the shape of the photon spectrum for wino searches. Here we provide the details of our NLL calculation, and perform a numerical study demonstrating that the theoretical uncertainty is significantly reduced when compared with the LL result, achieving an uncertainty from higher order corrections at the level of 5%. The cumulative cross section for z ≥ z cut is shown in Fig. 1 for a wino with a mass of 2.9 TeV, which corresponds to the case where the thermal relic density matches the measured one [32]. Here z cut restricts the cross section by allowing only photons with z ≥ z cut , see Eq. (2.2). We additionally show a band depicting the approximate values of z cut that correspond to the H.E.S.S. energy resolution. As can be seen, our calculation significantly reduces errors associated with the particle physics component of the annihilation cross section, which allows for the robust interpretation of experimental results in terms of DM model and astrophysical parameters. We also study the importance of using the full spectrum as compared with a single endpoint bin approximation when computing experimental limits by performing a mock H.E.S.S. analysis, and using the results of our forecasted H.E.S.S. limits [44]. We find that the use of the full spectrum near the endpoint is crucial to preserve the desired accuracy, emphasizing the importance of our EFT formalism.
Although the primary goal of this work is to provide a precision prediction for heavy WIMP annihilation in the endpoint region, a number of interesting features of SCET with broken gauge symmetry that have not previously appeared in the literature arise in our NLL calculation, and we devote several sections to their discussion.
An outline of this paper is as follows. In Sec. 2 we review the structure of the factorization formula for the endpoint region derived in [15], and prove that it remains valid at NLL accuracy. In Sec. 3 we discuss the necessary formalism for achieving resummation at NLL accuracy, give explicit results for all one-loop anomalous dimensions, and solve the relevant renormalization group (RG) equations. In Sec. 4 we present analytic results for the cumulative and differential spectra at NLL accuracy, and comment on some interesting aspects of their structure. Numerical results and a study of the related theoretical uncertainties are given in Sec. 5. In Sec. 6 we compare the use of the full spectrum with a single endpoint bin in a mock H.E.S.S. analysis and with our forecasted limits, emphasizing the importance of having a complete description of the shape of the photon spectrum in the endpoint region. We conclude in Sec. 7.

Factorization at LL order
In this section we provide a brief review of the factorization formula used to describe the photon spectrum in the endpoint region at LL order. A complete description including all notational conventions followed here, as well as a derivation of the formula through a multistage matching procedure, can be found in [15].
The essential process of interest is χ χ → γ X, where χ is the DM which annihilates to a hard photon γ, which is detected by the experiment, and additional radiation X. We will characterize the photon energy spectrum with the dimensionless variable We will be interested both in the differential spectra dσ/dz as a function of z, as well as the cumulative spectra as a function of z cut σ(z cut ) = 1 zcut dz dσ dz , (2.2) which is shown in Fig. 1. In the fully exclusive case (z = 1), only two bosons are produced in the final state. This is the relevant configuration for a line search with perfect resolution. However, given the finite energy resolution of real experiments, the region relevant for line searches is the so-called endpoint region, characterized by an energy resolution M χ (1−z cut ) M χ [15]. In this case, additional radiation beyond two bosons is present in the final state. To be near z = 1, this radiation must be either low energy, or collimated along the direction of the boson recoiling against the detected photon. This configuration is illustrated in Fig. 2a.
The collimated spray of radiation is referred to as a jet. Due to the phase space restrictions, large logarithms, log(M χ /m W ) and log(1 − z), appear in the perturbative calculation of the spectrum. In this paper we focus on improving the precision of the calculation relevant to the endpoint region. Our approach combines a number of different EFTs including Non-Relativistic DM and SCET as in [20][21][22], but with additional formalism to describe the photon endpoint energy spectrum. The formalism to describe the endpoint region developed in [15] makes use of the limit m W M χ (1 − z) M χ . This enables us to factorize the cross section into a number of different functions, each describing the dynamics at a particular scale. Log enhanced contributions to the cross section can then be resummed using RG techniques. The resulting resummed functions are then recombined to produce a precise prediction for the cross section.
As is well known, the Sommerfeld effect is relevant to heavy WIMP annihilation [27][28][29][30][31]. In our framework, the Sommerfeld effect can be factorized out of the cross section using where it is captured by the matrix elements In particular, we will require the following two matrix elements where χ 0 = χ 3 and χ ± = χ 1 ∓ iχ 2 / √ 2 . These matrix elements are obtained by numerically solving the associated Schrödinger problem using the approach detailed in App. A of [57]. We include a mass splitting, which for winos we take to be δ = M χ ± − M χ 0 164.4 MeV [58].
The focus of this paper is on extending the description of the radiation in the final state to a higher perturbative order. The starting point is the LL factorization theorem derived in [15]: where the dependence on the gauge indices a b ab is suppressed. Here we have explicitly added a superscript LL to emphasize that in [15] it was only shown that this factorization holds in this form at LL order. However, in the next section we will show that this formula holds also at NLL order, so that this higher order calculation can be achieved by improving the precision of the individual functions in Eq. (2.6). Here µ and ν are the virtuality and rapidity scales respectively, and we have defined the following shorthand notation for the convolution in the z variable (2.7) Operator definitions of the relevant functions will be given in Sec. 3 along with their anomalous dimensions. Here we restrict ourselves to providing a physical description of the different functions that appear in the factorization formula in Eq. (2.6). This factorization, along with the physical radiation described by each of the functions appearing in the factorization formula, is illustrated in Fig. 2b. The functions appearing in the factorization formula naturally divide up into two groups, those describing the dynamics above the electroweak symmetry breaking scale, which can be evaluated in the unbroken theory: • H M χ , µ : this hard function captures virtual corrections for χ χ → γ γ, γ Z.
and those that depend on the gauge boson masses and must be evaluated in the broken theory: • J γ m W , µ, ν : this photon jet function captures the virtual corrections to the outgoing γ.
• S m W , µ, ν : this soft function captures wide angle soft radiation at the electroweak scale.
• Jn m W , µ, ν : this jet function captures collinear radiation in X at the electroweak scale.
Overlap between the different functions is removed using the zero bin procedure [59].
In [15] each of these functions was computed to LL accuracy, and the consistency of the factorization was shown at that order. After showing that this factorization also holds at NLL, we will calculate each of these functions to NLL accuracy.

Validity of the Factorization at NLL
Various aspects of the LL factorization formula derived in [15] and displayed in Eq. (2.6) obviously remain valid at NLL (and higher) orders. This includes the factorization of hardcollinear and soft-collinear modes in the intermediate EFT, as well as the refactorization of the jet sector to separate m W M χ √ 1 − z. This suffices to define the functions H, J γ , Jn, and H Jn appearing in Eq. (2.6), but leaves a soft function S that is not fully factorized. Thus the key non-trivial aspect of the factorization formula in Eq. (2.6) is the refactorization of the soft function S , which describes low energy (soft) radiation.
In [15] it was shown that S could be factorized at LL into a hard matching coefficient H S , a collinear soft function C S (describing soft radiation collimated along the direction of the photon), and a soft function S. Explicitly, .
This refactorization is shown schematically in Fig. 3. It involves splitting the soft contributions into modes for S at the scale m W and modes for H S at the scale M χ (1 − z). These two modes are only separated in energy, but have no angular hierarchy, while in addition there are collinear-soft modes at the same energy scale as S but at a different angle. This causes a potential issue when proving Eq. (2.8), which no longer follows from a standard soft-collinear, hard-collinear, or hard-soft factorization. More precisely, soft emissions at the scale M χ (1 − z) are much more energetic than those at the scale m W , and therefore behave as eikonal sources for bosons at the scale m W . Any boson radiated at the scale M χ (1 − z) can radiate many bosons at the scale m W . An example of such a graph is shown by the red vertex in Fig. 4. From the perspective of the radiation at the scale m W , it then appears as though additional Wilson lines are present, beyond those in the original soft function, and one might be worried that more complicated Wilson line operators are required to capture all effects to NLL accuracy. This would imply that the simple factorized formula of Eq. (2.6) would need to be extended to go beyond LL order. It is known that this occurs in other cases, for example in the study of non-global logarithms (NGLs) [60] and their factorization [54,[61][62][63]. However, we will see that this behavior does not occur in the case studied here, and more complicated soft functions are not required at NLL.
We will show the validity of Eq. (2.8) at NLL through an explicit calculation. The refactorization in Eq. (2.8) asserts that one can describe the radiation by two separate soft functions, one describing radiation at the scale M χ (1 − z) and one describing radiation at the scale m W . A particular example of this factorization is shown in Fig. 4. To show that this is valid at NLL, one must show that there cannot be a logarithm arising from a graph where the bosons at the scale m W are color entangled with those at the scale M χ (1 − z). (In this  Figure 3. A review of the refactorization of the soft function derived at LL accuracy in [15]. In this paper we prove that this refactorization is also valid at NLL. (a) A schematic depiction of the refactorization of the soft function into different modes, illustrating their physical interpretation. (b) The virtuality and rapidity of the different modes appearing in the factorization formula. section we use the word "color" to refer to gauge index structure in the electroweak theory.) With two emissions, this corresponds physically to one emission at the scale M χ (1 − z) and one at the scale m W that did not factorize, as illustrated by the red vertex in Fig. 4. More precisely, this is the non-abelian component of the graphs, since eikonal emissions factor for the sum of abelian diagrams. Since the soft function is inclusive over radiation at the scale m W , one would naively expect that such logarithms do not exist, because of the cancellation between real and virtual corrections. However, due to the presence of electroweak charged initial and final state particles, the cancellation of IR divergences is not guaranteed by the Kinoshita-Lee-Nauenberg (KLN) theorem [64][65][66]. KLN violating effects, which arise due to a mismatch in color structures between the real and virtual corrections leading to an incomplete cancellation [67][68][69] could interfere with this argument. However, we can directly show that such logarithms do not exist by performing a two-loop calculation of the soft function in the strongly ordered limit, with one boson at the scale m W and one at the scale M χ (1 − z). The strongly ordered limit is sufficient to detect logarithms.
We now consider the calculation of the two loop soft function in the strongly ordered limit. We will separately consider two classes of diagrams: the triple gauge boson vertex and the independent emissions. We will show that for each class of diagrams contributing to a given soft operator, most of the real and virtual diagrams cancel pairwise. This is the case for all the triple gauge boson vertex diagrams as we discuss below. For the independent emission diagrams, the non-singlet nature of the soft function implies that the color structure of the real and the corresponding virtual diagram is not always the same. However, we still find that the non-abelian piece of the color structure cancels out, as is required to show the validity of To demonstrate the cancellation, we present explicit results for a particular soft operator, S 2 (whose precise definition is given in Eq. (3.41)), since it includes the most general structure involving all three types of Wilson lines in the n,n, and v directions. In our calculation, we will denote the two loop momenta by k 1 and k 2 , with k 1 ∼ M χ (1 − z)(1, 1, 1) being the harder and k 2 ∼ m W (1, 1, 1) the softer boson.

Triple Gauge Boson Vertex Diagrams:
We begin by considering diagrams with the triple gauge boson vertex. Representative real and virtual diagrams are shown in Fig. 5. These two diagrams are obtained simply by shifting the cut, and therefore their color structures are identical. In these diagrams, the higher energy particle has eikonalized from the perspective of the lower energy particle, and therefore this is precisely the type of diagram that one could worry generates an operator with an additional Wilson line structure. However, since we are fully inclusive over the low energy radiation, and the color structures of the real and virtual diagrams are identical, we will see that the real and virtual graphs cancel, and no logarithm is generated.
The cancellation in the strongly ordered limit is easily checked. Defining the real integrand as I R and the virtual integrand as I V , and taking the strongly ordered limit, we find Here F (k 1 ) is a function of k 1 whose precise form is not relevant for the current argument.
We can now perform the n · k 2 integral in I V by contours. Since k 1 crosses the cut, we have n · k 1 > 0. We can therefore close the contour in the upper half plane, and we only pick up the residue from 1/ k 2 2 − M 2 χ + i . This is exactly equivalent to the on-shell condition for k 2 with an overall minus sign. Thus the real and virtual contributions cancel. Using this same approach, it can be verified that this cancellation happens for all triple gauge boson vertex diagrams at two loops. Therefore, we see that for these graphs involving the triple gauge boson vertex, because the measurement function is inclusive over the low energy boson, no large logarithms are generated. We see explicitly that for these contributions, the presence of electroweak charged particles does not have an effect, since the color structure is manifestly identical between the real and virtual diagrams.

Independent Emission Diagrams:
We now consider the independent emission diagrams. Most of the independent emission diagrams cancel between the real and virtual diagrams that are related by a shifted cut, in exactly the same way as was just demonstrated for the triple gauge boson vertex diagrams. However, shifting in the cut does not always produce the same color structure for the real and virtual diagrams. We will decompose the color structure into an abelian and non-abelian piece, and we will then show that the non-abelian piece always cancels, as is required for our factorization.
We provide one illustrative example, which is shown in Fig. 6. Again, defining the real integral as I R and the virtual integral as I V , we find Here τ A and τ B are gauge factors from the k 1 and k 2 gauge boson-DM vertices respectively. To show that these two integrals cancel, we again perform a contour integral in n · k 2 . Closing the contour in the upper half plane, we only pick up the residue k 2 2 −M 2 χ +i , which implements the on-shell condition, similar to the triple gauge boson example. However, in this case the color structures are different, and therefore we do not have an exact cancellation. We can investigate this further by rewriting the color structure for I R to show that it is possible in both cases to separate these terms into an abelian and a nonabelian part, with identical color structure for the non-abelian contributions. Since the integrals have the same form with an opposite sign, these two expressions show that the nonabelian pieces cancel between the two diagrams, leaving behind an abelian term which describes the contribution to the soft function arising from independent emissions at the scale m W and M χ (1 − z). This contribution is already described by the factorization formula, which involves a soft function at each of these scales. This argument holds for all strongly ordered two-loop contributions to the soft function, allowing us to prove that it factorizes into two independent soft functions, one describing radiation at the scale m W , and one describing radiation at the scale M χ (1 − z). Therefore, in particular we have proven that the refactorized formula in Eq. (2.8) is also valid at NLL order.
Given the above argument, we find that our factorization formula for the resummed photon spectrum in the endpoint region can be extended to NLL without modification Again we suppress the gauge indices a b ab here. It would be interesting to understand if this factorization continues to hold at higher logarithmic orders. However, since it is not required for this paper, we do not pursue this further.

Renormalization Group Evolution at NLL
Having demonstrated the validity of our factorization formula Eq. (2.13) at NLL, deriving the cross section at this level of accuracy requires computing the anomalous dimensions of each of the functions appearing in Eq. (2.6) to O(α W ), and solving the associated renormalization group (RG) equations. After briefly reviewing the technology utilized to achieve NLL accuracy in Sec. 3.1, in Sec. 3.2 we provide results for all the one-loop anomalous dimensions appearing in our factorization, as well as the required matching coefficients. In Sec. 3.3 we then solve the associated RG equations. Due to operator mixing, and the large number of functions appearing in our factorization formula, this turns out to be non-trivial. This section is more technical than the remainder of the paper, and the reader interested only in the final result can skip it on a first reading.

RGE Preliminaries
We will resum logarithms appearing in the cross section using renormalization group equations (RGEs) in both virtuality, µ, and rapidity, ν. We will always choose an appropriate integral transformation of the factorization formula so that the RG is multiplicative. In particular, we will work in Laplace space, with s conjugate and As usual, the contour in the inverse Laplace transform is chosen such that γ is to the right of all singularities in the complex plane. The RGEs in virtuality will then take the form The anomalous dimension for the function F can be shown to have the following structure to the order required for our analysis Here Γ F andΓ F are proportional to the cusp anomalous dimension for the function F , which drives the double logarithmic evolution, and γ F [α W ] is the non-cusp anomalous dimension. To achieve NLL accuracy, the cusp anomalous dimension is needed to two loops, while the non-cusp anomalous dimension is needed to one loop. It is known that at two-loops that the cusp anomalous dimension Γ F [α W ] is a multiple of the universal cusp anomalous dimension Γ cusp [70], where c F andc F are constants depending on the function F , that do not have an expansion in α W , i.e., it can be determined using a one-loop calculation (and should not be confused with the fundamental Casimir C F ). We can expand the cusp anomalous dimension perturbatively as where here, and throughout the rest of the text we usẽ to simplify our notation. If the scale of α W is not made explicit, it should be taken to be evaluated at the natural scale of the function it appears in. The first two perturbative orders in the expansion of the cusp anomalous dimension, as required to achieve NLL accuracy, are well known [70] Γ 0 = 4 , (3.8) Here T F = 1/2 is the SU(N ) gauge group index, n f = 6 is the number of fermions in the fundamental representation, and C A is the Casimir for the adjoint representation. Note that here we have defined Γ cusp so as not to include an overall Casimir, which is included in c F andc F . We will similarly expand the non-cusp anomalous dimension inα W as We will also need to RG evolve in rapidity [71][72][73], and to do so we will use the rapidity RG framework of [72,73]. The rapidity RG equation takes the form where the anomalous dimension is given by As for the µ-anomalous dimension, Γ F ν (α W ) is a multiple of the cusp anomalous dimension, Γ cusp . Therefore, as with the µ-anomalous dimension, to achieve NLL accuracy, we need the one-loop non-cusp anomalous dimension, and the two-loop cusp anomalous dimension. We will find that the one-loop non-cusp contributions γ F ν vanish. At NLL, we will also need to take into account the running of the coupling, which is a single logarithmic effect. We define the perturbative expansion of the β function as where we will need where C F = 3/4 is the quadratic Casimir for the fundamental representation. The resummation of large logarithms is achieved by running all functions appearing in the factorization theorem to a single scale. Since the evolution in µ and ν commutes one can choose an arbitrary path in the (µ, ν) plane as long as large logarithms in the anomalous dimension are consistently avoided. In Fig. 7 we show the path that was used in [15] to simplify the RG evolution. Here all functions are run to the scale µ = m W . In principle, one must then run all the functions to a common ν value, which in the figure is shown as ν = 1/s. In practice, the ν running is trivial since the one-loop non-cusp rapidity anomalous dimension vanishes. Therefore, with this choice of path, it suffices to run the hard function and the hard matching coefficients for the jet and soft functions in µ, greatly simplifying the RG structure. Furthermore, we can treat the combination C S S =S as unfactorized, see Eq. (3.45) below. Beyond NLL, this would no longer be possible, which would significantly complicate the RG evolution.

Anomalous Dimensions and Matching Coeffecients
In this section, we give explicit results for all anomalous dimensions required to achieve NLL accuracy, along with the necessary matching coefficients. As discussed above, this requires the one-loop non-cusp anomalous dimensions. For brevity, we only give the final results, noting that all the required integrals can be found in App. A of [15]. Due to the large number of functions appearing in our factorization formula, instead of just giving the values of the cusp and non-cusp anomalous dimensions, we explicitly write the logarithms to show the natural scales appearing in the functions. Furthermore, to simplify the functions, we will write Γ cusp [α W ] with the assumption that when working to NLL, this cusp anomalous dimension should be kept to second order in the coupling.

Hard Function
We use a basis of hard scattering operators defined by [22] L (0) Here χ v are the non-relativistic heavy DM fields, which carry a label velocity v, as in heavy quark EFT [74,75]. We will take v = (1, 0, 0, 0). The B n⊥ are collinear gauge invariant fields [17,18], with labels n = (1, 0, 0, 1) andn = (1, 0, 0, −1). They are defined by Here D µ n⊥ is the collinear gauge covariant derivative, and W n is a collinear Wilson line where P µ is the label momentum operator, which when acting on a collinear field, returns its label momentum. The ultrasoft Wilson lines, Y , appearing in the operator are determined by the BPS field redefinition [18] B aµ n⊥ → Y ab n B bµ n⊥ , (3.18) which is performed in each collinear sector, and similarly for the non-relativistic DM particles. For a representation, r, we have where P denotes path ordering, and similarly for then direction. For our particular basis of operators, the relevant Wilson line structures are (3.20) At tree level, the Wilson coefficients are given by To NLL accuracy, the anomalous dimensions of the Wilson coefficients are given by [22] γ C (µ) = 2 γ W T (µ) 1 +γ S (µ) . Here the diagonal component contains the cusp anomalous dimension and the beta function, and there is also an off-diagonal non-cusp componentγ To simplify notation, we have defined which will appear frequently in our results.
The hard function is defined as We will use the simplified notation Note that a slightly different notation was used in [15], where we defined Beyond LL, it is necessary to treat H 12 and H 21 separately. For the same reason, when we redefine the soft operators below we will not follow the notation of [15]. 1 The hard function in this new basis satisfies an RG equation with anomalous dimension matrix which is given byγ The components of the anomalous dimension matrix are given by where where here and below we have explicitly used C A = 2 to simplify the formulas. At LL order this matrix is diagonal, while at NLL one encounters non-trivial operator mixing.

Photon Jet Function
The photon jet function J γ is defined as Its µand ν-anomalous dimensions are given by With our choice of resummation path, we do not need to RG evolve the photon jet function, and therefore only need its boundary value. However, to ensure a consistent definition of NLL accuracy in Laplace and cumulative space, it is well known that in Laplace space, one must also keep the O(α W ) RG generated logs in the boundary terms (see e.g. [76] for a detailed discussion). We will do this throughout the forthcoming sections without further comment. For the photon jet function, we have Here s W = sin θ W is the sine of the weak mixing angle and the canonical value for µ 0 γ is m W .

Recoiling Jet Function
The recoiling jet function is defined before refactorization by is the collinear measurement function, which returns the value of (1 − z) when acting on a collinear state. The recoiling jet function is refactorized into a low scale function Jn(m W , µ, ν) and a high scale matching coefficient The low scale function is defined as and is understood to be evaluated in the broken theory. Eq. (3.36) then defines the hard To NLL accuracy, the low scale function has neither a µor ν-anomalous dimension The hard matching coefficient only has a µ-anomalous dimension, which at NLL order is The boundary values for the matching coefficient and the jet function are given at tree level by

Soft Function
Before refactorization, we have the soft functions S where all functions S should be read as carrying gauge index structure a b ab. These definitions involve the measurement function which returns the value of (1 − z) when acting on a soft state. See [15] for more details.
Here, and in all subsequent expressions, we keep the time ordering and dependence on x = 0 implicit. Again note that the definition of the operators differs slightly from that used in our LL calculation [15], since we distinguish S 12 and S 21 at NLL. At NLL the non-vanishing µ-anomalous dimensions are given by The ν-anomalous dimension is diagonal in gauge index space, and is given bŷ We write the refactorization of the high scale soft function into a hard matching coefficient, a collinear-soft function, and a low scale soft function, as (3.45) As discussed above, we are able to choose a path such that we do not need to separately run the C S and S functions. Therefore, as in [15] we only give the anomalous dimension for the combined functionS, as well as for the matching coefficient H S . This significantly simplifies the structure. However, even with this simplification, we will find that there are 5 relevant refactorized soft functions, so that j runs from 1-5, and i runs from 1-4 in Eq. (3.45). This implies that there are 20 hard-soft functions. Fortunately, we will find that 10 of them vanish at this order, simplifying the structure of the RGE. For the low scale soft functions, we havẽ For simplicity, we have not written the free gauge indices on the left hand side of these equations. Here X n and V n are collinear soft Wilson lines and We have also defined the collinear soft measurement function, which returns the value of (1 − z) on a collinear soft state, see [15] for more details. This then defines the matching coefficients H S,ij . Again, we reiterate that as compared to [15] we have slightly modified our basis to incorporate the complete set of gauge index structures that are generated at NLL. The soft functionS satisfies a ν RGE where the matrix is diagonalγS The µ-anomalous dimension forS has a non-trivial mixing structure, The non-vanishing terms of this matrix are given by The matching coefficient for the soft function H S only has a µ RGE, which again exhibits non-trivial operator mixing. Making the indices explicit, we have Recall that here i runs from 1-4, and j runs from 1-5. We therefore have a system of 20 coupled differential equations. Fortunately, to NLL accuracy 10 equations drop out, due to the fact that   where ψ was defined in Eq. (3.25). These anomalous dimensions can be written in the form of a matrix equation as Due to the size of the matrices, we do not explicitly give them here, although they can be directly read off of Eq. (3.57). This provides the complete set of anomalous dimensions required to achieve NLL accuracy. We have checked that they satisfy all RG consistency constraints.

Solution to NLL Evolution Equations
Armed with all the necessary anomalous dimensions, we next turn to solving the RG equations at NLL. Due to our choice of resummation path, as was discussed in Sec. 3.1, we must µevolve all functions from their natural scale to the scale m W . Recall that the homogeneous µ-evolution equation takes the form where γ F µ (µ, ν) takes the form given in Eq. (3.4). For the purpose of the solution the terms not involving an explicit log(µ 2 ) can be grouped together, so for simplicity we use γ F [α W ] below, instead of writing out γ F [α W ] +Γ F [α W ] log(ν 2 /μ F (s) 2 ). Here we have suppressed all kinematic arguments for simplicity. At NLL we need to include the effects of the running coupling, so it is useful to change variables from µ to α W The solution to the RGE is then given by 2 where the evolution kernels are and To NLL accuracy, we have where recall Γ 0 and γ 0 control the cusp and non-cusp parts of the evolution and Using this evolution kernel, we can now run each of the functions to the scale m W . We will separately consider the jet function, hard function, and hard-soft matching coefficients. For the hard and hard-soft functions, due to the complicated mixing structure, we will have to first diagonalize the system of equations before using the kernels given in this section.

Recoiling Jet Function Evolution
The jet function must be evolved from its initial natural scale µ 0 J = 2 M χ √ 1 − z to the common scale m W . In Laplace space, the natural scale is which we have written in terms of where Θ is the Heaviside step function. The step function enforces that the mass of radiation in the jet is greater than m W . In addition to the RG evolution, we also need the appropriate initial value of the hard jet function at the initial scale µ 0 J . We remind the reader again that to ensure a consistent NLL accuracy in Laplace and cumulative space, in Laplace space one must also keep the O(α s ) RG generated logs in the boundary terms. For H Jn , we therefore have where C J contains the additional logs. Combining these results, we can then write down the hard-jet function evolved to the common scale m W as

Hard Function Evolution
We now consider the RG equation for the hard function. The hard function satisfies an evolution equation with non-trivial mixing at NLL where the explicit form of the entries were given in Eq. (3.29). We can write this as a diagonal evolution involving the Γ cusp and the beta function plus an off-diagonal non-cusp contribution Here Γ H is as defined in Eq. (3.30) and ψ is as defined in Eq. (3.25). The remaining non-cusp anomalous dimension can now be diagonalized where the matrix V defines the change of basis In this basis, writing Γ H explicity, we have the diagonal equations which we can now easily evolve from the natural scale µ 0 H (which has the canonical value 2 M χ ) down to the scale m W . Writing the non-cusp piece of the anomalous dimension as −4 β 0 + a, we can now define a hard evolution kernel where we have also defined Note that at the canonical scale the ω H contribution will vanish, but it is important to retain in order to estimate the impact of scale variation. We can then write Before inverting, we need the boundary values of the hard functions at their natural scale, which are given by where we used the results for the one loop cusp contribution to the hard scale matching coefficients from [25,82]. Substituting these in we find Note that although ψ appears in these results, the final result for the cross section will be purely real as it must -ψ will lead to the appearance of the sine or cosine of a phase. This occurs already for H 1 (m W ), and we see that a cosine of a phase appears in the resummed result. We will explain the physical origin of these phases in Sec. 4.1.

Soft Matching Coefficient Evolution
The soft matching coefficient H S satisfies an RG equation involving a 10 × 10 mixing matrix. To diagonalize the system, we perform the following invertible change of basis (3.82) After performing this change of basis, we obtain the following set of decoupled equations We can now solve these equations to evolve the functions from the high scale µ 0 S (with canonical value 2 M χ (1 − z)) down to m W . To do so, we define the soft evolution kernel which is given in terms of and all other boundary coefficients are zero, we obtain Here to simplify the expressions we have defined Recall that the remaining ten H S,ij (m W ) functions not listed here are zero. This solution resums all logarithms appearing in the hard-soft function to NLL accuracy.

Analytic Resummed Cross Section at NLL Accuracy
We now have all the relevant pieces to derive analytic cumulative and differential endpoint spectra for wino annihilation at NLL accuracy. The NLL cross section is given by where LP −1 denotes the inverse Laplace transform as defined in Eq. (3.2). In addition to previously defined shorthand notation, we have also defined the phases c H = cos 6 π β 0 log r H , s H = sin 6 π β 0 log r H , × C J and keeping only the terms relevant at NLL order. Note that Λ d = C J , where we have chosen to redefine it to make the notation consistent. In detail, we have where in analogy with L S (s) we have (4.5) The complexity of this result is both due to the multiple gauge index structures in the hard and hard-soft functions, and their contractions with the Sommerfeld factors.
To perform the inverse Laplace transform analytically, we set scales in cumulative space. The natural scales of the functions in Laplace space are therefore taken to be formally independent of the Laplace space variable s. The only required transform between Laplace space and cumulative space is . which appear in the Λ i (s) expressions, we use that in Laplace space all of these terms appear multiplied by an expression of the form s q . We can therefore rewrite these logs in terms of derivatives using ∂ n q s q = s q log n s .
The derivatives can then be evaluated analytically once the Laplace transform has been performed. These differential operators generate logarithms and polygamma functions. For detailed results see App. A. Using Eq. (4.1) and these steps we can then derive the following expression for the cumulative spectrum, as defined in Eq. (2.2): Here, for simplicity, we have only given the result evaluated with all scales at their canonical values. In addition we no longer show an s argument on the Λ i expressions, as all logarithms have now been replaced by derivatives according to Eq. (4.9). Results for general scales, as are required to generate scale variations, are given in App. A. The differential spectrum in z is generated by differentiating the result in Eq. (4.10) by z. We find: Here σ NLL exc is the NLL exclusive cross section, which is given by Our result for σ NLL exc reproduces the original NLL result for the wino from [22]. The analogous line result for scalar DM at NLL was computed in [21]. This provides a closed form NLL result which simultaneously resums all logarithms of (1 − z) and m W /M χ and correctly incorporates Sommerfeld enhancement effects. Eqs. (4.10) and (4.11) are the main results of this paper. Note that in Eq. (4.11) the Θ functions cut off the 1/(1 − z) power law singularities. Although this result is considerably more complicated than the corresponding LL expression presented in [15], it simply dresses the same 1/(1 − z) power law growth with additional logarithms. This structure persists to all logarithmic orders.

Non-Vanishing Electroweak Glauber Phase
Here we briefly comment on an interesting effect appearing in our final resummed result that occurs when electroweak charged external states are present. This effect is not specific to the case we are considering here, and indeed it also appears in fully exclusive calculations with charged electroweak states in the initial or final states (for the heavy dark matter annihilation case, this includes [21,22,82], and in the collider context from an EFT perspective [48-50, 83, 84]). However, recent advances in the understanding of the treatment of Glauber gauge bosons in SCET [85] give a clean interpretation of these terms. Since this connection has not previously been emphasized, we briefly deviate from our main goal to explain it here.
Our final result for the cross section involves the following phases, c H = cos 6 π β 0 log r H , s H = sin 6 π β 0 log r H . (4.13) These first appear at NLL in both the hard function H and the soft function S. These phases arise because the ψ = 1 − i π and ψ * = 1 + i π do not fully cancel in the cross section. This lack of cancellation has a physical origin when understood in terms of Glauber gauge boson contributions. In our factorization, the soft function describes soft modes with a homogeneous scaling n·p,n·p, p ⊥ soft ∼ Q λ 2 , λ 2 , λ 2 , (4.14) corresponding to on-shell degrees of freedom in the EFT. However, when evaluating the virtual soft integrals, we also integrate over regions of phase space where with a + b < 2. This is the scaling for the so-called Glauber region, which corresponds to off-shell exchanges that are instantaneous in both lightcone times. In our calculation, we do not treat the Glaubers as being separate from the soft modes. The fact that they are captured by the soft function in a hard scattering is due to their "cheshire" nature [85]. Once properly defined and regulated, Glaubers do not cross the cut, nor do they connect the initial to the final state. Example diagrams are shown in Fig. 8, illustrating the final state Glauber bursts G. These off-shell Glauber contributions to the virtual corrections give rise to the i π appearing in the one-loop amplitude, Eq. (4.13). These i π contributions at the level of the amplitude are ubiquitous. What is interesting here is that a phase contribution survives at the cross section level, since in many cases Glauber effects cancel for inclusive cross sections.
To understand why we are left with a phase in the cross section, we consider calculating the Glauber contribution to the S 12 soft function (the soft function arising from the interference of the two distinct gauge index structures). Computing the two relevant diagrams + , (4.16) which have Glaubers on either side of the cut, we find which is non-zero. Here the superscript G denotes the Glauber contribution. The two diagrams give different gauge index structures, which sum to the result shown here. Eq. (4.17) gives the i π terms that are included in our soft function S, which are related to the ones in H by renormalization group consistency. Note that if the external states were electroweak singlets, electroweak charge conservation would imply that the two diagrams would have the same gauge index structure. Then the two diagrams would be exactly conjugate, and the imaginary Glauber contribution vanishes, ψ + ψ * = Re[ψ]. In Eq. (4.17), this can be seen by contracting the result with a gauge index singlet final state δ ab S G,a b ab However, for wino annihilation the electroweak charges of the initial and final states are nonsinglet, and Glauber phases contribute to the cross section. Furthermore, the Glauber i π multiplies a logarithm, which once resummed yields the phases in Eq. (4.13). This can be viewed as a manifestation of the KLN violation in the Glauber sector, and is a completely general phenomenon when one has multiple electroweak charged initial or final states. This cancellation (or lack thereof) extends to multiple Glauber exchanges. It would be interesting to investigate the properties of the (non-) cancellation of electroweak Glaubers further. This is beyond the scope of this work since for our current application these phases appear in the hard-soft matching coefficient, and are correctly captured within our framework.   Figure 9. The cumulative (left column) and z 2 weighted differential (right column) spectrum for wino annihilation in the endpoint region. Spectra are shown for three different wino masses, 3, 10 and 35 TeV, at both LL and NLL accuracy. Theoretical uncertainties obtained by scale variation are shown by the shaded bands -the uncertainty bands are hardly visible at NLL at high masses.

Numerical Results and Uncertainties
In this section, we present numerical results using our NLL formula. In particular, we focus on the reduction in the theoretical uncertainty as compared with LL. The uncertainty bands are generated by scale variations, which probes higher order logarithms. Due to the structure of the µand ν-anomalous dimensions, as described in Sec. 3, we are able to choose a path up to NLL accuracy that does not require rapidity evolution (see Fig. 7). Since no logarithms at NLL are generated by ν-evolution, when performing the scale variation at LL, all logarithms at NLL can be probed by µ-scale variations. The LL uncertainty bands are generated by varying the µ-scales about their central value by a factor of 2.
The NLL error bands probe next-to-next-to-leading logarithmic (NNLL) logarithms. Note that at NNLL one has a non-trivial rapidity evolution between the collinear-soft and soft functions, so that they could not be combined into a single function as was done here. To capture this in our uncertainty estimate, our NLL uncertainty band requires both a µand ν-variation: explicitly we move both scales independently up and down by a factor of 2 for the photon jet function, and take the maximal band. This uncertainty is then added in quadrature with those from the other µ-scale variations. We believe that this is a reasonable estimate of the scale uncertainties. As a reference, we will find that the scale uncertainties for our result are comparable to those for the fully exclusive NLL wino calculation in [22], which also found 5% perturbative uncertainty. When the O(α W ) fixed order corrections were added in [25] to obtain NLL' accuracy, the perturbative uncertainties shrank to ∼ 1%, and the central value was near the boundary of the ∼ 5% NLL uncertainty bands. We expect the O(α W ) fixed order corrections to have a comparable effect in our calculation (some of these corrections, such as those to the hard function, are even identical). This further supports that our estimate of the higher order perturbative error is reasonable.
In Fig. 9 we show the cumulative (left column) and differential (right column) spectrum for wino annihilation for M χ = 3, 10, and 35 TeV (top, middle, and bottom rows respectively). This range of values was chosen since on the high end H.E.S.S. is currently probing DM masses up to 70 TeV [3], while on the low end our EFT expansion breaks down as we approach 1 TeV. For M χ 1 TeV, one should smoothly match our EFT onto an EFT where m W ∼ M χ (1 − z). This would be particularly interesting to consider for the Higgsino, which we leave to future work. We emphasize that the value 3 TeV (or more precisely 2.9 TeV) is particularly motivated as this is the mass that corresponds to a thermal relic wino. In all cases, we find that moving from LL to NLL yields a relatively small change in the central value. What is particularly noteworthy is that NLL demonstrates a large reduction in the theoretical uncertainty as compared to LL. With the NLL results in hand, the spectrum for heavy wino annihilation near the endpoint is now under excellent theoretical control.
The impact of this calculation is that it can be used to quantitatively explore constraints on the parameter space for winos from indirect detection. To this end, we provide Fig. 10, which shows the results of the mock analysis developed in [15], updated using our NLL prediction -more discussion on the details of the analysis is provided in the next section. For  [2], see [15] for details of the procedure. The comparison is made in terms of the line annihilation cross section σv line = σv γγ + 1 2 σv γZ (see text for details), as a function of the wino mass. Plotted here are the prediction (blue) and the mock limit (orange); the parameter space above the limit line is excluded assuming the Einasto DM profile. In both cases the bands represent the theoretical uncertainty associated with the NLL calculation. The overall normalization error is captured by the prediction band. For the mock limit, the uncertainty originates both from the variation in the shape of the endpoint spectrum and its normalization relative to the line. Finally, the thermal wino prediction M χ = 2.9 ± 0.1 TeV is also shown.
concreteness, the comparison between the mock exclusion and the theory prediction is made in terms of the line annihilation cross section σv line = σv γγ + 1 2 σv γZ . We note that in order to convert from the endpoint cross section σ NLL computed here to σ line , one must evaluate Eq. (4.11) in the limit z → 1, and be careful to keep track of the fact that here we are computing the rate for γ + X, which introduces a factor of 2 in the conversion since both of the γ's from σ γγ contribute, i.e., and the approximation made here treats the kinematics for γ γ and γ Z identically, see [15,44] for additional discussion of this convention. Note that these mock limits only include the contribution from the line and endpoint spectrum; the justification to neglect the contribution from continuum production resulting from wino annihilation to W + W − at lower masses was provided in [15]. While we caution that a genuine analysis of the 2013 H.E.S.S. data should be done to provide an actual limit, we see that our mock limit shows that the thermal wino with /s]

Thermal Wino Cross Section
Endpoint Inclusive Semi-Inclusive Exclusive Figure 11. A comparison of our LL and NLL calculations with inclusive, exclusive and semi-inclusive predictions from the literature. The best agreement is found with the semi-inclusive calculation. The disagreement as z cut → 1 is due to unresummed logarithms of (1−z cut ) in the semi-inclusive calculation which are correctly captured using our formalism.
mass M χ = 2.9 TeV is excluded by a factor of ∼ 25. We also emphasize that this assumes an Einasto DM profile, so that one way to avoid this seemingly stringent bound on wino annihilation is to core the profile, see e.g. [44,57] for a discussion. Importantly, the theory error band shown for the prediction in this plot is now under excellent control, justifying the need for our NLL calculation. For a more careful exploration of the implications of the NLL endpoint spectrum in the context of H.E.S.S. forecast limits, see [44]. Finally, in Fig. 11 we show a comparison of our NLL cross section with several calculations that exist in the literature. In particular, we compare with the fully exclusive (line) calculation at NLL [22,25], the inclusive calculation at LL [20,23], and the semi-inclusive calculation at LL [24]. With the reduced NLL uncertainties, we see that for z cut ∼ 0.8-0.9, our prediction differs significantly from the exclusive and inclusive predictions, being approximately intermediate between the two, which individually each sum large log(M χ /m W ) logarithms at NLL order. As expected, the semi-inclusive provides a better approximation, agreeing with the shape and norm of the LL endpoint result away from z cut → 1. However, this calculation does not resum logarithms of 1 − z cut , which become important as z cut → 1. The effect of these logarithms, which are properly captured in our result, is clearly seen by the fact that the semi-inclusive result rapidly diverges from our LL and NLL endpoint results as z cut → 1.

Endpoint Spectrum Versus Fixed Bin Approach for Experiments
Obtaining a reliable theoretical interpretation of indirect detection line searches requires correctly incorporating experimental constraints into the underlying theoretical setup. This consideration has resulted in the appearance of a number of different approaches: fully inclusive [20,23], fully exclusive [21,22,25], semi-inclusive [24], non-zero fixed bin width [82], and endpoint spectrum [15]. Clearly, what is observed by the experiment is the true photon spectrum convolved with the appropriate instrument response functions. In our endpoint approach, this can be correctly treated since we have computed the full shape of the photon spectrum. With the NLL calculation presented here, we have a theoretically robust result with an estimated 5% residual perturbative uncertainty.
In this section we investigate how well the full result could be approximated by assuming that the resolution effects at an experiment such as H.E.S.S. are captured by integrating the photon spectrum over a bin with some effective width (with the bin maximum being z = 1). This is important since it allows for an understanding of which approximations are valid for correctly describing the experimental setup. One motivation is that in certain cases determining the rate in a single bin near the endpoint is easier than achieving the full spectral shape. If this approach were demonstrated to be a good approximation, it could also pave the way for more calculational efficiency.
To address this question we compare the experimental constraints on the wino obtained using our full spectrum to those derived by rescaling constraints on a gamma-ray line. One might hope that for an appropriate choice of bin size, the true constraint would trace the line limits up to a simple rescaling. However, we will demonstrate that there does not appear to be any simple way to choose such a bin a priori.
In the interest of thoroughness, we will provide the results for two different analyses (for more details, see [15] and [44] respectively) 1. Mock Analysis: A simplified re-analysis of a 2013 H.E.S.S. result [2], following the methodology we developed previously [15]. We perform a χ 2 analysis on the data, floating a seven-parameter background model in addition to the signal. The functional form of the background model was chosen [2] to provide a good description of the data in the region of interest. We confirm that we approximately reproduce the quoted constraints on a pure line, and then apply the same analysis taking our full NLL spectrum as the signal. We estimate the energy resolution for this dataset based on interpolating the values given for photon energies of 0.5 and 10 TeV [2].
2. Forecast Analysis: A detailed forecast of H.E.S.S. sensitivity, based on Monte Carlo simulations of the expected background, the instrument response functions, and the analysis pipeline [44]. We perform a likelihood analysis that incorporates both spatial and spectral information, for both signal and background. In this section we show results (reproduced from [44]) assuming an Einasto profile for the DM density, although this method is generalizable to alternate DM density profiles. We compute the sensitivity  We see that in both approaches, the bin width required to map the line-only signal onto the limit derived for the full spectral shape is a non-trivial function of the DM mass.
to a pure line and compare that to the prediction of our NLL spectrum. This forecast includes cuts on the data that allow an effective energy resolution of 10% independent of energy.
We begin by understanding the effect the various approximations would have on H.E.S.S. limits, in Fig. 12a (Fig. 12b) we plot the mock (forecast) limits for σv line , considering a number of different ways of modeling the experimental analysis given the theory calculation. First, in solid orange we show the limit derived using the NLL prediction, which incorporates the full spectral shape. The black dotted line shows the limit assuming the pure line analysis using the resummed exclusive prediction for the annihilation rate. As was shown in [15], the inclusion of the photon distribution near the endpoint almost always enhances the limits. Next we consider a number of different approximations for effective bin widths. In dashed green, we show what happens if instead of using the shape of the distribution, we take our NLL prediction for the spectrum and just assume a bin of width equal to the H.E.S.S. resolution (as appropriate for the calculation in question). Although this is not a terrible approximation, it is well outside our estimate of the theoretical uncertainty at NLL -this emphasizes the need for a correct treatment of the spectrum when working at this accuracy. Finally, we also consider the limits derived when taking a bin width of m W /M χ and (m W /M χ ) 2 . These scales are natural from the EFT perspective, as they tie the bin width to the other physical scales in the problem, simplifying the calculation [20,23,82]. The appropriate EFT for a bin width m W /M χ uses the scaling m W ∼ M χ (1 − z), while the appropriate EFT for a bin width (m W /M χ ) 2 uses the scaling m W ∼ M χ √ 1 − z. Relative to our EFT, these include additional power corrections. However, at the level of accuracy that we work, namely NLL, we do not include the matrix elements in the low energy theory, only logarithms from evolution, and therefore we should correctly reproduce the large logarithms in these bins for the m W ∼ M χ (1 − z) case, which also bounds the m W ∼ M χ √ 1 − z case. Therefore, for simplicity, we do not perform dedicated analyses for these bins, and instead simply take our EFT result evaluated with bin widths of these values. We believe that the qualitative conclusions that we draw should be robust. Results for a resolution of (m W /M χ ) 2 are shown by the dashed blue curves, and for realistic H.E.S.S. resolutions we see that this bin size is too small (with our approach the dashed blue curves fully overlap the line curves, but will presumably differ somewhat in a full treatment of m W ∼ M χ √ 1 − z as in [82]). Finally we see that a EFT setup which utilizes a m W /M χ bin width gives the dashed purple curves, which lie in-between the line and endpoint results for M χ 3 TeV. These EFTs are perhaps more relevant for a DM mass of around 1 TeV, e.g. in the case of the Higgsino. However, the region of DM masses they can describe is likely to be rather limited since they are forced to have a resolution that scales like 1/M χ . This is distinct from the behavior of the H.E.S.S. resolution, which is approximately flat with increasing DM mass.
Another aspect of Fig. 12a and Fig. 12b that deserves comment is the fact that the mock limit and the forecast limit performances are similar for the important mass range near M χ ∼ 3 TeV, even though the forecast limit assumes a larger data set. This behavior occurs for two reasons. First, the forecast limit utilizes a more sophisticated background determination method, with less flexibility to incorrectly absorb a signal into the background model, and less sensitivity to fluctuations in the background determination. Consequently, estimating the relative strength of the expected limits in the two analyses is non-trivial. Additionally, the mock analysis is based on observed data, not expected data, and the 2013 dataset from which the mock limit was derived placed a stronger-than-expected upper limit on the flux near 3 TeV [3], likely due to a downward fluctuation in the background.
Another way of testing the fixed-bin approach is to explore the extent to which the H.E.S.S. endpoint analysis can be reproduced by using a single endpoint bin with an effective size, and compare that size with the H.E.S.S. resolution. To this end, we provide results computed by rescaling the limit on (or sensitivity estimate for) the line cross section, proportionally to the number of photons in the endpoint bin, and then determine what bin size is needed to reproduce the constraint (or sensitivity estimate) obtained by the endpoint analysis involving the full spectrum. In Fig. 12c, we plot the effective bin size required to reproduce the constraint obtained by the endpoint analysis in the first calculation above, compared with the true H.E.S.S. resolution for the 2013 dataset. We see that the effective bin size varies non-trivially as a function of the DM mass, which could not be predicted without a full calculation of the shape of the distribution. Indeed, at low masses the limit with the full spectrum is weaker than the limit from the line only. This is likely due to the non-trivial interplay between the signal and the flexible background model used in the mock analysis (based on that of Ref. [2]). These results emphasize the need to convolve the NLL spectrum computed here with the experimental line shape and full background fit. To further demonstrate this point, Fig. 12d shows the same result for the sensitivity estimate from the second calculation described above, using stateof-the-art background modeling methods. The requisite effective bin size is somewhat more stable with mass in this case; sensitivity estimates do not include the statistical fluctuations inherent to real data, and the background modeling also has less freedom in this analysis. However, the bin size needed to match the predicted level of sensitivity again does not match the nominal H.E.S.S. energy resolution over the vast majority of DM masses, and shows non-trivial variation with DM mass.
One might also ask whether further improvements in the accuracy of the endpoint spectrum would allow corresponding improvements in the experimental constraints. For the H.E.S.S. experiment, however, the theoretical uncertainties at NLL are currently subdominant to the experimental systematic uncertainties. In more detail, the Galactic Center (GC) region is a very crowded environment in very-high-energy (VHE) gamma rays. The determination of the residual background, mostly coming from misidentified CR hadrons, is challenging in the GC due to the presence of gamma-ray sources and regions of diffuse gamma-ray emission. A proven robust approach to probe signals from a cusped DM density distribution in the GC region consists of masking regions of the sky with known VHE emission, and then making use of the reflected background method where background and signal are measured in the same observational and instrumental conditions allowing for a precise determination of the background level [3,86]. However, systematic uncertainties arise from the imperfect knowledge of the energy scale and energy resolution of the instrument. In addition, the Night Sky Background (NSB) rate -corresponding to the unavoidable optical photon light emitted by bright stars in the field of view -varies significantly over degree scales in the GC region. This variation induces a systematic uncertainty in the background determination when using the reflected background method. Propagating this uncertainty into the DM constraints implies a systematic uncertainty in the limit on the line annihilation cross section ranging from a few percent up to 60%, depending on the DM mass [3], which dominates over the O(5%)

NLL theoretical uncertainties.
Future studies will make use of precise Monte Carlo simulations for the expected residual background determination in the GC, in the same observational and instrumental conditions as for the signal measurement. This could alleviate the level of systematic uncertainties substantially; for example, the inhomogeneous NSB could be accurately simulated in each pixel, allowing a careful subtraction of this component. The main remaining systematic uncertainty might then become the uncertainties in the energy scale and energy resolution of H.E.S.S.; a systematic uncertainty in the energy scale of 10% shifts the limits by up to 15%. The energy resolution uncertainties are a smaller effect; while the energy resolution is weakly dependent on the observational conditions, and this can impact the limit, a deterioration of a factor of two in the energy resolution only induces a decrease of 25% in the expected limit [3]. Greater precision in the theoretical calculation thus might eventually become valuable from an experimental perspective, but would require multiple current sources of systematic uncertainty at the 10% and higher level to be reduced below 5% by improved analyses.
The results of this analysis demonstrate that to have a reliable theoretical interpretation of the experimental results requires a calculation of the full shape of the photon spectrum. Approximations using effective bin widths are simply not reliable at the level of accuracy of an NLL calculation. In the event that multiple sources of experimental systematic error can be reduced in the future, extending the calculation presented here to NLL+NLO accuracy (or NNLL) could become necessary.

Conclusions
In this paper, we have extended the calculation of the hard photon spectrum for wino annihilation in the endpoint region, as is relevant for indirect detection experiments, to NLL accuracy. This calculation was performed using an EFT framework developed in [15], which facilitates the factorization of distinct physical effects. In particular, our result includes both the resummation of Sudakov logarithms and the Sommerfeld enhancement. The theoretical uncertainties of our calculation are of the order of 5%, which is a significant reduction as compared with our earlier LL prediction. In particular, the theory uncertainties are now sufficiently under control so as to make a subdominant contribution to the total uncertainty relevant for experimental exploration.
In the course of our calculation we encountered a number of interesting effects associated with electroweak radiation and the presence of electroweak charged initial and final states. For example, we found that the non-electroweak-gauge-singlet nature of the incoming and outgoing states led to a non-trivial remaining Glauber phase in the final cross section. We think this would be interesting to explore further in a more general context than that considered here.
The H.E.S.S. telescope has collected a large dataset of photons from the Galactic Center region, with an energy resolution of ∼ 10%, permitting sensitive searches for spectral features. Using both a mock H.E.S.S. analysis and a detailed forecasting framework, we studied the importance of using the full photon spectrum computed in this paper as compared with a single endpoint bin approximation when computing experimental limits. We find that the mapping to an effective bin width is a highly non-trivial function of the DM mass. This emphasizes the importance of having theoretical control over the shape of the distribution in the endpoint region, and not simply the photon count in an endpoint bin, for deriving accurate limits from experimental data.
With an understanding of our factorization formula at NLL accuracy, it is now straightforward to calculate the spectrum at this accuracy for other heavy WIMP candidates, such as the pure Higgsino, the mixed bino-wino-Higgsino, or the minimal DM quintuplet [29,[87][88][89][90][91]. This will allow for the robust theoretical interpretation of indirect detection constraints for these compelling DM candidates from the wealth of data at current and future experiments.