Energetic $\gamma$-rays from TeV scale dark matter annihilation resummed

The annihilation cross section of TeV scale dark matter particles $\chi^0$ with electroweak charges into photons is affected by large quantum corrections due to Sudakov logarithms and the Sommerfeld effect. We calculate the semi-inclusive photon energy spectrum in $\chi^0\chi^0\to \gamma+X$ in the vicinity of the maximal photon energy $E_\gamma = m_\chi$ with NLL' accuracy in an all-order summation of the electroweak perturbative expansion adopting the pure wino model. This results in the most precise theoretical prediction of the annihilation rate for $\gamma$-ray telescopes with photon energy resolution of parametric order $m_W^2/m_\chi$ for photons with TeV energies.


Introduction
Within the large variety of dark matter (DM) candidates a weakly interacting particle with mass m χ in the 100 GeV to 10 TeV range (WIMP) stands out due to its conceptual minimality and its relation to the electroweak scale. Although loop suppressed, the pair annihilation of WIMPs into two photons or a photon and a Z boson often provides a distinctive signature in the form of a monochromatic component of high-energy cosmic γ-rays. The strongest limit on such γ-line signals are currently set by the H.E.S.S. experiment [1][2][3]. This limit is expected to be improved by an order of magnitude by the Cherenkov Telescope Array (CTA) [4] under construction, which will severely constrain generic WIMP models. This motivates a thorough theoretical investigation of the basic pair annihilation cross section given a particular model.
It is well-known that TeV-scale DM annihilation is not accurately described by the leading-order annihilation rate, but is modified by the Sommerfeld effect [5] generated by the electroweak Yukawa force on the DM particles prior to their annihilation. In terms of Feynman diagrams the Sommerfeld effect corresponds to ladder diagrams with W , Z, photon and Higgs boson exchange, which contribute O((m χ α 2 /m W ) n ) at order n, where m W is the W boson mass and α 2 the SU(2) gauge coupling. For annihilation rates to exclusive final states, in addition to the Sommerfeld effect large logarithmically enhanced quantum corrections of O((α 2 ln 2 (m χ /m W )) n ), also known as electroweak Sudakov logarithms, arise [6] from the restriction on the emission of soft radiation. 1 Electroweak Sudakov logarithms in DM annihilation into photons have been identified as a potential source of large corrections and resummed to all orders in perturbation theory in previous work [7][8][9][10][11] in models with an electroweak triplet scalar or fermionic DM particle.
In this letter we revisit this question starting from the observation that telescopes do not measure two photons from a single decay in coincidence. Rather, the observable will be an enhancement of photons in a single bin, whose energy corresponds to the maximal photon energy E γ = m χ of the semi-inclusive photon spectrum γ + X, where X denotes the unidentified other final-state particles of the annihilation process. While the leading term in the expansion in α 2 is indeed exclusively from the γγ and γZ final state, the logarithmically enhanced quantum corrections are different for the exclusive and the semi-inclusive measurement. If E γ res m χ denotes the energy resolution of the instrument for photon energies E γ ≈ m χ , the maximal invariant mass of the unobserved final state X is m 2 X = 4m χ E γ res m 2 χ . Since X must balance the large momentum of the observed photon, X is forced to have a 'jet-like' structure. The kinematic situation involving non-relativistic heavy particles in the initial state and energetic, small-invariant mass objects in the final state is naturally described by a combination of non-relativistic and soft-collinear effective field theory, similar to the QCD treatment of the 'inverse' situation of hadronic production of two heavy particles [12].
In the following sections we first summarize briefly the basic elements of an effec-tive field theory (EFT) treatment of the single-inclusive photon spectrum d(σv)/dE γ in DM pair annihilation near the kinematic endpoint. For the DM model we refer to the widely discussed pure wino model, which features an electroweak triplet whose electrically neutral component is the DM particle, although some results apply more generally to DM particles in an isospin-j multiplet. We then present and discuss our result for the all-order resummed spectrum including both the Sommerfeld and Sudakov corrections. A more detailed exposition of the formalism as well as extensions will be reported in a longer article. While this work was being finalized, a similar EFT calculation of the endpoint of the γ + X spectrum has appeared [13]. The present EFT formulation refers to a finer photon energy resolution, but includes the one-loop corrections to all matching coefficients, soft and jet functions thus achieving NLL' rather than NLL accuracy for the observable in question.

The resummed energy spectrum
We add to the Standard Model (SM) Lagrangian a fermionic multiplet χ (which can be of Majorana or Dirac type) in an arbitrary isospin-j representation of the electroweak (EW) SU(2) gauge group. For the Majorana case, only integer j are allowed, while for the Dirac case also half-integer j are possible. In both cases we assume zero hypercharge (Y = 0). The DM particle is the electrically neutral member χ 0 of the 2j +1 dimensional multiplet. The Lagrangian is when χ is a Dirac fermion. For the Majorana case, χ is self-conjugate and its Lagrangian is multiplied by 1/2. The SU(2) covariant derivative is D µ = ∂ µ − ig 2 A C µ T C where T C , C = 1, 2, 3, are the SU(2) generators in the isospin-j representation and A C µ are the EW gauge bosons. In these models the dark matter particle obtains the correct relic density from thermal freeze-out for m χ in the 1-10 TeV range [14] for the favoured small representations j = 1 2 , 1, 3 2 , 2.

Effective theory framework
We consider the process for nearly maximal photon energy. Since the kinetic energy of the dark matter particles is negligible, E γ max = m χ . Assuming an energy resolution E γ res of the γ-telescope, we are interested in the quantity which together with the astrophysical line-of-sight factor determines the flux of photons from dark matter annihilation into the energy bin that contains the photon line signal. The constant c ∈ [0, 1] accounts for the fact that m χ may not coincide with the upper energy value of the bin. The photon endpoint spectrum depends on four important scales: m χ (hard), the small invariant mass m X = 4m χ E γ res (collinear) of the unobserved, energetic final state, enforced by the kinematics of the endpoint, the electroweak scale m W (soft) and the energy resolution scale E γ res (ultrasoft). We shall now assume that the energy resolution is parametrically of order E γ res ∼ m 2 W /m χ , which implies m X ∼ m W and the scale hierarchy E γ res m W , m X m χ . The factorization of the multi-scale Feynman diagrams into single-scale contributions, which is a prerequisite to all-order resummation, then requires the introduction of momentum modes with the following parametric scaling: Here λ = E γ res /m χ and k µ ∼ (n + · k, n − · k, k ⊥ ) where n µ + , n µ − are two light-like vectors with p µ γ = E γ n µ + and n + · n − = 2. We remark that the collinear, anti-collinear, soft and potential modes all have the same virtuality O(m 2 W ). The interactions of these modes are described by standard potential non-relativistic and soft-collinear effective Lagrangians (similar to [12], generalized from QCD to the electroweak interaction).
One might also consider the wider resolution E γ res ∼ m W , which implies E γ res , m W m X m χ and a different mode structure. Conceptually, the main difference is caused by the fact that the previous, narrower resolution does not allow the radiation of soft particles with electroweak-scale masses into the unobserved final state. Although the resolution of the up-coming γ-ray telescopes is probably closer to the wide resolution case, in this work we concentrate on the narrow resolution E γ res ∼ m 2 W /m χ to stay close to the line-like signal. The wide resolution case, which is in fact simpler from the EFT point of view, will be discussed in subsequent work, which will also provide the explicit forms of the effective Lagrangians.

Factorization
The primary annihilation process is described at leading order in an expansion in λ, which is also an expansion in m W /m χ , by operators O i for the S-wave annihilation of the dark matter particles into two EW gauge bosons. Once the hard modes are integrated out into the coefficient functions C i of these operators, the collinear, anti-collinear and potential fields can no longer interact directly, since their momenta would add up to hard virtualities. The collinear modes build up the unobserved final state X, while the anti-collinear modes must result in the single, observed photon. The non-relativistic DM particles are described by the potential fields exchanging potential EW gauge bosons, which causes the Sommerfeld effect.
The factorization formula then follows from an analysis of the coupling of the (ultra) soft modes to any of the above. The decoupling of soft gauge boson attachments from heavy-particle fields in the presence of the Sommerfeld effect by a time-like Wilson-line field redefinition of the heavy-particle field has been demonstrated in [12] for an unbroken gauge theory. Although in the present case the Sommerfeld effect must be computed in the broken theory with gauge boson masses, soft attachments still factorize from the ladder diagrams, since a soft momentum throws the potential heavy particle off-shell, which removes the enhancement of the ladder rungs between the soft attachment and the hard vertex. It follows that the Sommerfeld factor S IJ completely factorizes from the Sudakov resummed annihilation rate Γ IJ (E γ ), where the sums over I, J run over all electrically neutral two-particle states that can be formed from the 2j + 1 single-particle states of the electroweak DM multiplet. For example, in the triplet ('wino') model, I, J = χ 0 χ 0 , χ + χ − . Since gauge boson exchange between the DM particles prior to annihilation can change the initial two-particle state χ 0 χ 0 into any I, J, the annihilation rate with the Sommerfeld effect factored out is a matrix describing the amplitude for the annihilation of two-particle state I times the complex conjugate of the annihilation amplitude for state J. The Sommerfeld factor is defined in terms of the matrix element of non-relativistic DM fields of the schematic form χ 0 χ 0 |[χχ] J |0 0|[χχ] I |χ 0 χ 0 such that S IJ = δ IJ in the absence of the potential force causing the Sommerfeld effect (see [15] for the full definitions). The coupling of soft gauge fields to collinear and anti-collinear fields is removed by a redefinition of the (anti) collinear fields with light-like Wilson lines [16]. Since the small energy resolution forbids soft radiation into the final state X, the soft function is a vacuum matrix element of Wilson lines that can be regarded as a soft Wilson coefficient D of the annihilation amplitude [17,18] on top of the hard Wilson coefficients C i of the operators O i . At this point all modes have been factored into separate functions except for the coupling of ultrasoft fields. At leading order in the expansion in λ, the coupling of ultrasoft gauge fields is again described by Wilson lines, where now the gauge field must be the photon. Since the initial state is electrically neutral and at rest, the coupling of ultrasoft photons is of higher-order, leaving ultrasoft Wilson lines from the (anti) collinear fields. There is no kinematic restriction on the radiation of ultrasoft particles with energy and masses of order E res γ into the final state, hence the ultrasoft function corresponds to an amplitude squared summed over the unobserved ultrasoft final state.
We can therefore write down the following factorization formula for the energy spectrum of DM annihilation (more precisely, the off-diagonal annihilation matrix Γ IJ in the DM two-particle states) into a single identified-photon inclusive final state near the maximal photon energy: The prefactors account for the spin average, flux factor and photon momentum angular integration, C and D refer to the hard and soft matching coefficients of the annihilation amplitude. The former is evolved from the scale 2m χ to the electroweak scale µ W . The indices V, W, X, Y refer to the SU(2) adjoint representation of the electroweak gauge bosons. In the second line Z γ is the anti-collinear factor for the observed photon, J the jet function for the unobserved collinear final state X, convoluted with the ultrasoft function S. Since the (anti) collinear and soft modes have parametrically equal virtualities, a rapidity evolution factor V (µ W , ν s , ν j ) is needed. We postpone the derivation of this formula and technical details to a separate paper and proceed by defining the appearing functions and providing the results of their calculation.

Operator basis and hard matching coefficients
The relevant hard annihilation operators can be written as Here χ v is a two-component non-relativistic spinor field in the SU(2) weak isospin j representation, χ c v = −iσ 2 χ * v the charge-conjugated field, and A B ⊥c,µ the collinear gaugeinvariant collinear gauge field of soft-collinear effective theory (SCET). Collinear fields have large momentum component n + p. A similar definition applies to the anti-collinear direction with n + and n − interchanged. Fields without position arguments are evaluated at x = 0. The operators are non-local along the light-cone of the collinear directions. T A are SU(2) generators in the isospin-j representation. The spin matrix is defined by (conventions n µ ± = (1, 0, 0, ∓1), 0123 = −1) The first form holds in d dimensions in dimensional regularization, but it turns out that evanescent operators are not important, since the non-relativistic, soft and collinear dynamics is spin-independent. n is a unit vector in the three-direction. Since the (anti) collinear gauge field in the operator is transverse, the Lorentz indices µ, ν and corresponding spatial indices m, n are effectively always transverse. Note that the DM bilinear is in a spin-singlet configuration.
We normalize the effective annihilation Lagrangian as The one-loop calculation of the MS-subtracted matching coefficients results in whereĝ 2 (µ h ) is the SU(2) gauge coupling in the MS scheme at the matching scale µ h ∼ 2m χ , and c 2 (j) = j(j + 1) the SU(2) Casimir of the isospin-j representation. 2 The coefficients are evolved to the EW scale µ W µ h . The evolution is diagonal [19] in the basis O i where the DM bilinear transforms in an irreducible SU(2) representation given by such that The evolution factor in the irreducible isospin-J representation satisfies the renormalization group equation with anomalous dimensions 2 For the wino model j = 1, the one-loop coefficients were previously given in analytic form in [11] in the context of resumming the annihilation rate to the exclusive γγ, γZ final state. When transforming to their operator basis, we find that our coefficient ofĝ 4 2 /(16π 2 ) differs by +4 (−4) from their C 1 (C 2 ). We could track this difference to an error in the external field renormalization for the DM field and an inconsistency in combining counterterms for Dirac and Majorana χ fields. We note that the final result for the coefficients functions is independent of whether the DM particle is a Majorana or Dirac fermion.
γ J H,s (α 2 ) = γ H,s c 2 (J) Here c 2 (ad) = 2 is the Casimir for the adjoint representation and n G = 3 denotes the number of fermion generations. (The Higgs contribution −16/9 to the two-loop cusp anomalous dimension has been obtained from the -scalar contribution in [20].) The NLL' and NLL approximations require the SU(2) cusp anomalous dimension in the SM in the two-loop approximation, and the other anomalous dimensions at the one-loop order, as given explicitly above. The LL approximation makes use only of the one-loop cusp term and neglects the other anomalous dimensions. In Figure 1 we show the evolved coefficient functions in the above mentioned approximations. The evolution equation is solved by numerical solution of the differential equation in the given approximation after solving the coupled system of renormalization group equations for the three gauge couplings, the top Yukawa and the Higgs self-coupling in the two-loop approximation.

Soft functions
The soft renormalization factor of the annihilation vertex is given by the vacuum expectation value of the Wilson lines that arise from decoupling soft EW gauge bosons from the fields that appear in the operators O i . In (6) the soft factor D i I,V W is defined in the basis of DM two-particle states with respect to the SU(2) indices of the DM bilinear, which corresponds to the definition The matrix K takes the linear combination appropriate to the two-particle state I, and arise from the gauge fields and are in the adjoint representation. The time-like Wilson line Y v is in the isospin-j representation of the DM field. All four Wilson lines extend from x = 0 to infinity in their respective directions. Note that although V, W = 1, 2, 3 are the gauge boson indices referring to W A rather than the mass eigenstates W ± , Z, γ, the soft function lives at the electroweak scale and must be computed with the Feynman rules of the SM after electroweak symmetry breaking, including gauge boson masses, contrary to the hard coefficient functions discussed above, which can be computed in the unbroken theory, neglecting the masses of the SM particles.
The requirements of an observed energetic photon and electric charge conservation imply that only the index values V, W, X, Y = 3 contribute to (6), so the corresponding sums disappear. The NLL' approximation requires the one-loop calculation of every function that appears in the factorization formula. For the triplet ('wino') model we find for the two operators and the two distinct two-particle states I = 00, +−. Since there are three regions with equal virtuality O(m 2 W ) but different light-cone momentum component n + · k, the soft function is not well defined with only a dimensional regulator. We use the rapidity regulator [21] in addition to dimensional regularization to obtain the above result and the jet functions below. The scale ν is related to the rapidity regulator. The soft function contains no large logarithms if µ, ν ∼ O(m W ).

Photon jet function
The 'jet' function for the exclusive anti-collinear photon state is defined by the squared matrix element of the three-component of the transverse SU(2) gauge field. A B ⊥µ denotes the gauge field dressed with anti-collinear Wilson lines, can be interpreted as the on-shell photon field renormalization constant in light-cone gauge.

Jet function of the unobserved final state
The jet function pertaining to the inclusive (unobserved) collinear final state is defined as the total discontinuity of the gauge-boson two-point function Again the field A B ⊥µ refers to the collinear gauge-invariant gauge field, which equals the ordinary gauge field in light-cone gauge.
We compute the 33 component in Feynman gauge to the one-loop order and write J 33 (p 2 ) in the form The one-loop correction is split into two contributions, of which the first refers to diagrams that involve at least one contraction with a gauge field from a Wilson line in the definition of A B ⊥µ and the second to the remaining diagrams, which are of self-energy type, as shown in the first and second line of Fig. 2, respectively. Only the Wilson-line diagrams require the rapidity regulator, and their sum is given by where The self-energy type contribution J 33 se (p 2 ) can be expressed in terms of conventional oneloop gauge-boson self-energies, which can be found, for example in [22]. We have taken the massless-quark limit of these expressions except for the top quark. We will provide the lengthy expressions in the detailed write-up.
In Fig. 3 we show the dependence of the integrated jet function on the invariant mass of the unobserved collinear final state. The integrated function jumps from a value aroundŝ 2 W (m W ) to a value around 1 as the invariant mass passes through m Z , and then slowly increases. The range shown contains the W + W − , ZH and tt thresholds (m t = 173.2 GeV is used here), which, however, are barely visible. The singularity of the NLO correction near m Z can be removed by a proper treatment of Z boson resonance. However, below we will adopt values E γ res > m 2 Z /(4m χ ), which implies that p 2 max is always larger than m 2 Z . The jet functions contain no large logarithms when µ = O(m W ) and ν = O(m χ ). The 'rapidity logarithms' related to the different value of ν that minimizes the logarithms in the soft and jet functions, respectively, can be summed at NLL' by solving the rapidity renormalization group equations [21] in the one-loop approximation. For the case at hand we find We checked that in the sum of all contributions the poles in the dimensional and rapidity regulator cancel. The hard, soft and jet functions above are defined by minimally subtracting the poles.

Ultrasoft function
The kinematics of the process does not allow soft radiation with momentum of order m W into the final state, which prohibits EW gauge boson radiation. However, radiation of photons and light quarks with masses of order or less than m 2 W /m χ is possible, which implies the convolution of the unobserved-final state jet function with an ultrasoft function accounting for the energy taken away from the collinear final state by ultrasoft radiation. The ultrasoft function is defined in terms of Wilson lines of ultrasoft photons and depends on the electric charges and directions of the particles in the initial and final state. After factoring the Sommerfeld effect, also the χ + χ − initial state must be considered. But for the S-wave annihilation operators O i only the total charge of the initial state is relevant for ultrasoft radiation, which vanishes. Furthermore, only the electrically neutral 33 components of the (anti) collinear functions appear for the γ + X final state. We therefore conclude that the ultrasoft function is trivial, S γ (ω) = δ(ω) . For this reason did not indicate the ultrasoft scale dependence of the functions in (6), and the convolution integral in (6) disappears.

Sommerfeld factor
The various functions discussed above are assembled according to (6) into the annihilation matrix Γ IJ of the χ 0 χ 0 and χ + χ − DM two-particle states. The photon energy spectrum is then obtained according to (5) by tracing this matrix with the Sommerfeld factor S IJ . For the fermionic DM triplet model, the Sommerfeld factors were first computed in [5]. In the present work we employ the modified variable phase method [15] for solving the Schrödinger equation and, different from the above, use the on-shell value α 2 = 0.0347935 of the SU(2) coupling as well as α OS (m Z ) in the Yukawa-Coulomb potential for the two-state system. Near the Sommerfeld resonance, the result is very sensitive to the mass splitting between the electrically charged and neutral members of the triplet. The mass splitting with two-loop accuracy can be inferred from [23,24], and is given after adjusment to our input coupling parameters by The first and second Sommerfeld resonances are located at 2.285 TeV and 8.817 TeV, respectively, for the coupling parameters and mass splitting employed in this work.

Results
It is straightforward to calculate the one-dimensional integral (3) that defines the γ + X yield from DM pair annihilation in a photon energy bin of size E γ res . For the discussion below we shall assume 4m χ E γ res = (300 GeV) 2 , which implies that the unobserved final state includes γ, Z, W + W − , ZH and light fermion pairs in the collinear jet function. We also adopt c = 1 in (3), since the resolution in an actual experiment will almost certainly be worse than assumed here. Figure 4 shows (upper panel) our results for σv (E γ res ) as defined in (3). The displayed DM mass range includes the first two Sommerfeld resonances. The four lines refer to the Sommerfeld-only calculation, which employs the tree-level approximation to Γ IJ (black-dotted) and the successive LL (magenta-dashed), NLL (blue-dashed) and NLL' (red-solid) resummed expressions for the same quantities with the latter representing the best approximation.
The importance of resummation of electroweak Sudakov logarithms becomes more apparent by normalizing to the Sommerfeld-only result (lower panel). As expected and already seen in previous LL and NLL calculations [7][8][9][10][11] of related exclusive and semiinclusive final states, resummation reduces the annihilation rate and increasingly so at larger DM mass. In the interesting mass range around 3 TeV where wino DM accounts for the observed relic density, the rate is suppressed by about a factor of two.
The resummed predictions are shown with theoretical uncertainty bands computed from the separate variations of the scales µ h , ν j in the interval [m χ , 4m χ ] and the scales µ W , ν s in [m Z /2, 2m Z ], added in quadrature. We find that at the NLL' order the residual theoretical uncertainty from scale dependence is negligible -in the figure it is given by the width of the red line. Our result can be compared most directly with the recent work by Baumgart et al. [13] who considered the same γ + X final state at larger resolution E γ res m W with NLL accuracy. We observe that the inclusion of one-loop corrections to the hard, soft and jet functions in our NLL' computation has the main effect of eliminating the theoretical uncertainty of resummation by reducing the scale dependence from 23% (LL) to 3% (NLL) to 0.4% at NLL' at m χ = 2 TeV. A similar reduction of scale dependence was already observed in the NLL' calculation of the exclusive γγ, γZ final state [11]. In numbers we find that at m χ = 2 TeV (10 TeV), the ratio of the resummed to the Sommerfeld-only rate is 0.536 +0.128 −0.121 (0.285 +0.104 −0.085 ) at LL, 0.579 +0.032 −0.003 (0.322 +0.018 −0.002 ) at NLL and 0.572 +0.004 −0.000 (0.318 +0.002 −0.000 ) at NLL'. The central values are evaluated at the central scales of the intervals above. We also varied all four scales simultaneously within these intervals and determined the maximal variation. The scale dependence at NLL' with this more conservative procedure increases by about a factor of two, which does not change the general picture.
In conclusion, we computed the γ + X spectrum near maximal photon energy from electroweak triplet ('wino') DM annihilation including the resummation of the Sommerfeld effect and electroweak Sudakov logarithms in the NLL' order. The inclusion of the electroweak one-loop corrections at NLL' renders the theoretical uncertainty of resummation negligible. It is plausible that the dominant theoretical uncertainty now arises from the fact that the non-relativistic EFT is only employed at leading order in the computation of the Sommerfeld enhancement, and from O(m W /m χ ) power corrections. In a subsequent paper we shall present an extension of this work to the case of wider photon resolution, further details on the EFT framework, the derivation of the factorization formula, and a comparison with expected experimental limits. The computations performed here are presently restricted to simple DM models, which add to the SM a single electroweak multiplet. It would be of interest to extend them to more complex models such as the MSSM, which would put the analysis of indirect detection constraints for mixed DM models [25] on the same theoretical footing as for minimal models.
Physics' (SFB 1258) and the Excellence Cluster 'Universe' of the Deutsche Forschungsgemeinschaft. MB thanks the Albert Einstein Institute at Bern University and the Department of Energy's Institute for Nuclear Theory at the University of Washington for hospitality during the preparation of this work.