Constraints on Galactic Wino Densities from Gamma Ray Lines

We systematically compute the annihilation rate for neutral winos into the final state gamma + X, including all leading radiative corrections. This includes both the Sommerfeld enhancement (in the decoupling limit for the Higgsino) and the resummation of the leading electroweak double logarithms. Adopting an analysis of the HESS experiment, we place constraints on the mass as a function of the wino fraction of the dark matter and the shape of the dark matter profile. We also determine how much coring is needed in the dark matter halo to make the wino a viable candidate as a function of its mass. Additionally, as part of our effective field theory formalism, we show that in the pure-Standard Model sector of our theory, emissions of soft Higgses are power-suppressed and that collinear Higgs emission does not contribute to leading double logs.


I. INTRODUCTION
One of the biggest challenges in contemporary high energy physics is determining the identity of dark matter (DM). Unfortunately, all of the astrophysical and cosmological ev- The WIMP Miracle presents a particularly compelling link between the weak scale and dark matter (see [1]). Demanding the correct relic abundance from cosmological freeze-out leads one to an O(TeV)-mass particle with electroweak-strength coupling. We consider here the case of one particular DM candidate, the wino, that belongs to a supersymmetric explanation of the weak scale. Despite the lack of direct evidence for the MSSM, the discovery of a SM-like Higgs at 125 GeV leaves open the possibility of a modestly-tuned supersymmetry scenario that retains a simple mechanism of SUSY-breaking and a standard, thermal-relic DM particle [2]. In fact, the only feature of the MSSM we use is the presence of a stable, electroweak triplet fermion. Thus, our result does not depend on the larger supersymmetric story, but holds for any DM scenario with the same quantum numbers, annihilating primarily through its gauge interactions. Extensions to scalar triplet or other SU(2) representations are straightforward.
A nearly pure wino offers one of the simplest supersymmetric DM candidates. In theories of anomaly-mediated SUSY-breaking, it emerges as the lightest supersymmetric particle (LSP) [3]. Furthermore, thermal-relic bino dark matter generically overcloses the universe by several orders of magnitude. If we assume that the wino constitutes all of the dark matter -which is an assumption we will relax below -and that its relic density was set at freeze-out, then the mass is constrained to the window M Wino ≡ (M χ ) = 2.7-2.9 TeV [4,5].
In principle, we may further constrain the wino (DM) via direct detection. However, the cross section for a TeV wino to scatter off nucleons is σ ∼ 10 −46−48 cm 2 , putting it far below current limits [6]. However, since the wino can annihilate directly to photons, by searching for monochromatic, O(TeV) photon lines, we can hope to discover it via indirect detection.
The authors of [4,7] used limits from the HESS Cherenkov telescope to argue that the nonobservation of such a photon feature put wino DM in severe tension with experiment. 1 In particular [7] calculated the annihilation rate to be ∼ 15× larger at M χ = 3 TeV than the HESS limit.
However, as with any indirect detection experiment, we must take into account astrophysical uncertainties. In particular, [4,7] consider variations of the DM halo profile for the galaxy. While cuspy profiles are preferred by simulations, it is possible that the dark matter density flattens out to a "core" about the galactic center, and such a distribution would lead to fewer DM annihilation events along our line of sight to the galactic center. To alleviate the tension with HESS, some amount of coring will be necessary to return pure wino dark matter to viability, but the question, which we answer in Section VI, is how much?
The annihilation rate of two neutral nonrelativistic particles, cannot be reliably calculated at tree level since it is plagued by infrared (IR) divergences which are cut-off by the gauge boson mass, M W ∼ 100 GeV. These divergences manifest themselves as large radiative corrections of two types. 2 One set comes from the potential interactions of the slowlymoving DM and scales as powers of α W Mχ M W > ∼ 1; the resummation of these corrections results in a Sommerfeld enhancement to the rate [8]. This effect can increase the rate by as much as O(10 4 ) relative to a perturbative calculation and is therefore a crucial step in analyzing wino DM. The second type of IR sensitivity is a Sudakov double-log, α W log( 2 , that can enter inclusive observables due to the non-singlet nature of our external states. This effect is known as "Bloch-Nordsieck Theorem Violation" and is generically found in the Higgs phase of non-Abelian gauge theories [9,10], such as the electroweak sector. Computations of fixed, NLO corrections to the exclusive, two-body annihilation rate found a 75% reduction relative to tree level plus Sommerfeld enhancement [7,11]. This opened up the possibility that the wino could still be viable, even with a non-cored dark matter profile. This result motivated the need for a systematic approach to the calculation of the rate since such a large radiative correction gives the appearance of series which is diverging too soon for an asymptotic expansion with such a small coupling. In a previous paper, we derived a factorization theorem that was used to derive the leading-log (LL) semi-inclusive wino annihilation rate (χ 0 χ 0 → γ + X) [12] in terms of a 1 They also consider the constraints from continuum gamma ray emission provided by the Fermi experiment.
These are most useful for constraining low-mass, few-hundred GeV winos. As the radiative corrections we investigate are much weaker in this regime, we do not investigate it here. 2 The term "divergences" is used despite the fact that the rate is physical. model-dependent set of matrix elements. The semi-inclusive rate is the relevant observable for constraining the wino with HESS since only a single hard photon from annihilation is measured and the resolution of the experiment (∼ 400 GeV at 3 TeV) is too poor to distinguish two-body from n-body annihilation (e.g. χ 0 χ 0 → γ + W + W − ) [13]. Despite the fact that both the Sommerfeld and Sudakov effects arise from the same hierarchy, M χ M W , they can be factorized through a mode decomposition of the relevant fields. Physically, we are separating the regime of the slowly-moving WIMPs evolving in the presence of the electroweak potential from the highly-energetic gauge bosons in the final state. We therefore employ a theory that combines a nonrelativistic treatment of the WIMPs, similar to NRQCD [14], with Soft-Collinear Effective Theory (SCET) [16] for the light annihilation products.
Hybrid theories of this form have appeared in the analysis of the photon spectrum in radiative decays of quarkonium [17] and electroweak SCET was elucidated in [18]. More recently, other groups performed the Sudakov resummation for WIMP annihilation, albeit for the exclusive two-body rate, employing the same effective field theory (EFT) [23,24].
In this paper, we give more details of of the factorization theorem presented in [12] in Section II and elucidate the role of the soft and collinear Higgs fields at higher orders in Section III. In Section IV we review the anomalous dimension calculations that allow Sudakov log resummation and provide additional detail on our use of the rapidity renormalization group [25]. To obtain a full calculation of the annihilation rate, in Section V we present our analysis of the Sommerfeld enhancement for the particular case of wino dark matter with a parametrically large Higgsino mass. We find that compared to tree level plus Sommerfeld corrected rate, the leading-log radiative corrections lead to a modest, few-percent reduction in the semi-inclusive rate to an observable photon. In Section VI, we also present exclusion plots for wino dark matter as a function of the mass and the amount of coring in the dark matter profile and the wino fraction of the dark matter. After concluding, we detail in an appendix the derivation of the quantum mechanical potential used in our Sommerfeld enhancement calculation from the underlying, relativistic quantum field theory.

II. FACTORIZATION
To develop a factorization theorem it helps to work in an EFT, which in our case is a hybrid of SCET and NRQCD that power counts in a double expansion in v, the relative velocity of the DM particles and λ = M W /M χ . We will work to leading order in both these parameters, which seems quite reasonable until a discovery is made.
The nonrelativistic (NR) piece of the Hilbert space is described by an effective theory which is analogous to NRQCD, but differs in two important aspects. First off, our theory of interest is in the weakly coupled Higgs phase and moreover, the potential between the DM particles is screened by the gauge boson mass. This latter distinction changes the way in which we power count as we shall see below. For various power counting schemes in NR theories see [19]. Despite these distinctions, the modal analysis of the effective theory follows from the general discussion in [15] and furthermore we will still use the acronym NRQCD.
NRQCD can be formulated in terms of three distinct types of modes, each with a unique scaling of momenta: Potentials (E ∼ mv 2 , p ∼ mv), soft (E ∼ mv, p ∼ mv) and ultra-soft where v is the relative velocity of the massive nonrelativistic states. The massive states have energy-momenta scaling as potential modes but nonetheless are on-shell, due to their NR dispersion relation. On the other hand, the potential gauge bosons are off-shell and can be integrated out to form non-local potentials. The soft modes are necessary to generate the correct running of the potentials, but will not play a role in the theory at hand as all of their effects will be of higher order. The same can be said for the US modes. However, the US modes will play a crucial role in determining the gauge invariant structures allowed in the theory.
It is important to understand that v is NOT necessarily the incoming relative velocity of the NR particles. v can also be the virialized velocity. That is, if v is sufficiently small then the particles can inspiral, gaining kinetic energy until the system virializes such that V (r) ∼ mv 2 . In the Coulomb phase (or for sufficiently large source masses in the confining phase of an asymptotically free theory) we have the relation which leads to the scaling v ∼ α given the Bohr radius r ∼ 1/(αm). If the incoming relative velocity is large enough then the system may not virialize. The condition for virialization is that the leading order potential be non-perturbative. In the Coulombic case this condition is α v ∼ 1. When this condition is met we must sum the box graphs. In such cases the incoming relative velocity becomes irrelevant for the power counting. 3 The NRQCD-like theory in our case is more complicated since gauge boson exchange flips a neutralino to a chargino which is taken to be a few hundred MeV heavier. Moreover the charge state admits Coulomb exchange, although the off-shell nature of the chargino intermediate plays the role of an IR cut-off. Since we are only interested in working at leading order in the v, the exact details of the correct effective theory will be irrelevant. All that matters for the present analysis is that there exist some virialized velocity, v, which will play the role of the power counting parameter. Given that additional rungs on ladder diagrams bring inverse powers of v, we need to sum over an infinite ladder of potential mode gauge boson exchanges between the massive fermions resulting in the so-called Sommerfeld enhancement which we will discuss in Section V.
where we use the vacuum insertion approximation in the WIMP sector, which is valid up to O(v 2 ) corrections. Henceforth, we drop the explicit vacuum projector. Implicitly, there is also a projection onto a single-photon state between the B µ⊥ fields i.e.
where X contains the accompanying particles in the collinear jet. All operators which arise in the matching can be reduced to one of these four using the Majorana condition. The only relevant nonrelativistic bilinear isχγ 5 χ. The spin one operators are irrelevant since Fermi statistics would lead to an antisymmetric SU(2) initial state, and we are interested in the annihilation of two neutral particles. Furthermore, P-wave annihilation is velocity suppressed. We have also used the definition where the ⊥ symbol implies the component perpendicular to the large light cone momentum n · p, where n µ = (1, 0, 0, 1) and D ⊥ µ is the covariant derivative in the collinear sector (for details see [16]). This field interpolates for a collinear gauge boson and is invariant under collinear gauge transformations due to the placement of the two collinear adjoint Wilson lines defined by The collinear Wilson lines are also implicit in the fermionic χ fields, where they also guarantee manifest collinear gauge invariance, where ζ B n is the fermion obtained by expanding QCD in the n-collinear limit and integrating out the small spinor components.
For a general kinematic configuration, the χ fields do not transform under SCET-II soft gauge transformations, since such a transformation would throw the χ field off-shell. Thus naively it appears that these operators are trivially soft gauge invariant. Indeed this is true, but it does not mean that the soft mode does not play a role in this process. The soft contribution to the operators arises after one integrates out the off-shell intermediate states which arise when softs couple to collinear. The question then becomes, where should the soft Wilson lines be inserted into our operators? We could perform a matching calculation, however there is a much simpler method which we call the "method of descent", developed in [26]. Here we will use a variation of those arguments.
The idea is to choose a kinematic scenario in which the the invariant masses of the external states are such that the soft momenta in SCET-II and the US modes in the NRQCD sector have the same scaling. We also raise the virtuality of the collinear modes so that soft radiation leaves the collinear lines on-shell 5 . In this scenario soft gauge invariance uniquely 5 When there is a hierarchy between the invariant masses of the soft and collinear momenta the theory is usually called SCET-I and the soft fields are called ultra-soft. For the sake of clarity we will call all non-collinear fields in SCET "soft". fixes the positions of the Wilson lines. The invariant mass of the external states is then lowered to its physical value, keeping the soft Wilson lines fixed.
To apply this methodology to the case at hand we tune M χ such that M χ v 2 ∼ M W , and we raise the virtuality of the collinear modes to be of order p 2 c ∼ M χ M W . This would be appropriate if we were to, say, measure the jet mass and not the photon energy, in which case the color structure of the operator basis would remain as in Eq. 2. Physically, in this limit, US can communicate between the NR matter fields and the collinear modes of SCET because such interactions leave both modes on shell. This apparent breakdown of factorization is remedied by performing a BPS [16] field redefinition where there are two types of path ordered adjoint soft Wilson lines S v and S n defined by This field redefinition decouples the US and collinear fields at the level of the action and dresses the operators such that O 2 and O 4 become The operators O 1 and O 3 receive no soft corrections. We now continuously deform M χ back to its physical value. In doing so, the soft fields retain their invariant mass of order M W , and the soft Wilson lines remain fixed by continuity.
The annihilation spectrum may be written as where O a s = S T vA A S vBB S T nÃA S nBB (11) and F γ AB is a fragmentation function defined by and F γ = F γ AB δÃB. Note that this is an unusual fragmentation function in that we are measuring states which are not gauge singlets. However, electroweak symmetry breaking makes this physically meaningful since it is trivial for observers to agree on the appropriate, EM-preserving gauge.
C 1-4 are the matching coefficients that give the probability for the dark matter to annihilate and create a photon with momentum n · p. F γ is the canonical fragmentation function giving the probability of an initial photon with momentum k to yield a photon with momentum fraction n · k/n · p after splitting. Since the contribution in Eq. (10) proportional to F γ is not sensitive to the nonsinglet nature of the initial state, it will only contribute large double logs from mixing with O 2,4 .
In writing down Eq. 8, we factorized the collinear and soft fields, as the total Hilbert space of the system is a tensor product of the soft and collinear sector. In general, none of the NRQCD modes can interact with the SCET mode without throwing them off-shell, thus leading to power suppressed interactions. Of course, the size of the power corrections will be dictated by this offshellness, but independently of the system's details, these interactions will not lead to large double logs.

III. THE ROLE OF THE HIGGS
Given that the Higgs mass is of order the weak scale, in principle we should consider both soft and collinear Higgs emissions in our factorization theorem. However, as we shall now show, the couplings of both collinear and soft Higgses are power suppressed. Assuming that the Higgsino is sufficiently heavier than the winos we may neglect the couplings of the Higgs to the nonrelativistic sector of the Lagrangian. In cases with light Higgsinos, the potential is affected, as discussed in [31].
The coupling of the Higgs to the gauge bosons is given by Let us consider the coupling of a soft Higgs to collinear gauge bosons. Recall that for a collinear gauge field with large momentum n · p, the polarizations scale as where the power counting parameter is λ ∼ M W /M χ . The soft higgs field scales as λ, while the measure scales as λ −3 . 6 Thus the coupling of a soft Higgs to a collinear jet is down by one factor of λ.
Now suppose we are interested in the coupling of a collinear Higgs to a collinear gauge boson. In principle this could, as in the case of soft Higgs emission, lead to non-analyticity that must be reproduced by the effective theory. The coupling of collinear Higgs with collinear gauge bosons (in the same light cone direction) is leading order since the measure will now scale as λ −4 . Such interactions will be written in the effective theory as This interaction will generate running in the fragmentation functions at subleading orders only.
It is interesting to ask whether or not the emission of a collinear Higgs in the n direction and will be neglected.

IV. CALCULATING THE ANOMALOUS DIMENSION
Much of the analysis in this section we first presented in [12]. To calculate the anomalous dimensions we first introduce an operator basis in the collinear and soft sectors where there is an implicit sum over the polarizations of the photon. The operator O b s has a trivial structure and hence does not receive radiative corrections, meaning its anomalous dimension is 0. Since it is a color singlet operator, the real and virtual poles cancel in the corrections to O b c , thus its anomalous dimension is 0 as well. At tree level, the vacuum matrix element of the operator O a s is simply δ A Ã δ B B . At one loop, the diagrams that contribute to this matrix element are shown in Fig. 2. The soft and collinear modes have the same virtuality and hence the divergences that arise from the factorization of the soft sector from the collinear cannot be regulated by dimensional regularization, which respects boost symmetry. Hence, we need to introduce a rapidity regulator, which manifestly breaks boosts [25]. This requires a corresponding factorization scale which we call ν. Using this formalism for the soft sector gives us . (17) Thus, we see that even though O a s itself has a nontrivial color structure, at one loop it generates a color singlet piece in addition to nonsinglet. We can do a similar calculation for the vacuum matrix elements in the collinear sector (cf. Fig. 3).
These operators clearly mix within their respective sectors and we can define anomalous dimension matrices for the scales µ and ν.
The anomalous dimensions are thus given by γ c ν,aa = 3g 2 4π 2 log( Since the rapidity regulator was needed to handle the divergence that arose from our artificial cut between the soft and collinear sector, any trace of it, including dependence on the scale ν should vanish when we combine soft and collinear results. Therefore, the fact that the ν Notice that the RHS of Eq. 23 is independent of the rapidity scale as it must be.
To present a model independent form we have used the fact that the tree level result for χ 0 χ 0 → γ + γ/Z must vanish in which case there is only one independent matching coefficient. This gives us that C 1 = C 4 , C 3 = 0 and C 2 = −2C 1 for matching at the high scale M χ . Using these boundary conditions we can solve for the Wilson coefficents at the low scale µ ∼ M W .
Similarly for C 3 and C 4 By running the Wilson coefficients, we resum all the leading double log contributions.
The cross section can now be obtained by evaluating the effective theory matrix elements at their natural scale µ ∼ M W . Here, all matrix elements are merely their values at tree level.
At the low scale, we are working in the broken theory, where the mass eigenstates are the neutralino χ 0 and the charginos, χ ± , which are defined as Anticipating the form of the Wilson coefficient at the scale M χ , we pull out a delta function δ(E γ − M χ ). Thus, utilizing the simplifications above along with the new basis we obtain the final form of the cross section up to corrections in the relative velocity [12].
where f ± ≡ 1 ± exp[− 3α W π log 2 ( M W Eγ )], which we plot in Figs. 4 and 5. Because f + > f − by a factor of a few throughout our range of interest and the numerical prefactor and Sommerfeld factors, as we will discuss in the Sections V, maintain this modest hierarchy, one can obtain a quick estimate of the effects of radiative corrections and resummation to the rate by looking at Fig. 5. We define the wavefunctions as where |χ 0 χ 0 S = 1 √ 2 (|χ 0 ↑ (p 1 )χ 0 ↓ (p 2 ) − |χ 0 ↓ (p 1 )χ 0 ↑ (p 2 ) ), and here after we drop the "S" subscript since our nearly-static, annihilating Majorana fermions are automatically in the spin singlet. Writing in terms of two component fields makes the singlet structure manifest. This gives usχ 0 γ 5 χ 0 = 2(χ 0 ) T iσ 2 χ 0 where the two-component fermions on the RHS are the large-component subset of the fourcomponent spinors on the LHS in the nonrelativistic limit. In order to fix the Wilson coefficient C 1 , we match onto the tree level annihilation cross section of a spin singlet chargino state 1 ). The leading order cross section to γ + X (χ + χ − → γγ + 1 2 (χ + χ − → γZ)) is given as From the effective theory description, we calculate this particular cross section in the nota-tion of Eq. 28 as At tree level, f + = 2 and 0|( which gives us Comparing the two results we fix

V. SOMMERFELD ENHANCEMENT
In order to quantify the semi-inclusive rate calculation, we need to determine the wavefunction-at-the-origin factors that enter our final, LL-resummed differential cross section in Eq. 28. The wavefunctions themselves are defined in Eq. 29 and can be computed in principle in the nonrelativistic effective theory by summing the ladder exchange of electroweak gauge bosons between winos to all orders. Fortunately, this is equivalent to the operationally simpler task of solving the Schrödinger equation for our two, two-body states |χ 0 χ 0 and |χ + χ − in the presence of the electroweak potential [8,27,28]. In Appendix A, we detail the process of obtaining this potential from the underlying quantum field theory.
Since it contains Coulomb, Yukawa, and mass-shift pieces and is off-diagonal for the two states, we solve it numerically. As expected for slowly moving particles in the presence of an attractive potential, we find Sommerfeld enhancement for the annihilation, that for some regions of M χ is orders of magnitude above the perturbative rate.
Taking into account appropriate state normalization, the Schrödinger potential is where δM ≡ M χ + − M χ 0 . For numerical analysis, we use δM = 0.17 GeV, which is its value over much of MSSM parameter space. To solve the system, we also need appropriate boundary conditions for our wavefunctions ψ 1 ≡ r|χ 0 χ 0 (ψ 2 ≡ r|χ + χ − ), where r is the relative distance between the two particles in the state, and we always work in the center of mass frame. The total wavefunction for our two state system is ψ = (ψ 1 ψ 2 ). We are ultimately interested in the annihilation of the neutral state, which is controlled by physics at length scales ∼ 1 Mχ 1 M W , and is therefore quantified by the wavefunction-at-the-origin, ψ 1,2 (0). We thus see that for the boundary condition of an incoming neutral, spin-singlet state, in the notation of Eqs. 28 and 29, where the numerical factors account for the difference between our definition in Eq. 29 and the interpolating fields for the two-body states and the tree-level normalization of the matrix element below Eq. 32. The latter can be found by comparing the potential we use for Schrödinger evolution (Eq. 35) with the nonrelativistic field theory interactions we obtain in Eq. A5. Since the annihilation process is perturbative, we can calculate ψ 1,2 (0) by turning off the annihilation and solving the scattering problem for the electroweak potential. 7 This is equivalent to our factorized effective field theory setup (cf. Eq. 28), where the wavefunctions are computed in the nonrelativistic effective theory, but the annihilation process is given by the perturbative, high-energy Wilson coefficient.
To proceed, we adopt the general analysis of [30,31], and stick to the former's notation as much as possible, for our state of interest. As is generic for a scattering problem in the presence of a central potential, for our wavefunction, where we recognize the incoming plane wave and scattered, outgoing radial wave. The index n labels the charged or neutral component of the two-body state, with n = 1 being neutral.
Since our asymptotic, physical state is a pair of neutral winos with CM frame E = Mχv 2 4 , we demand that the incoming plane wave only be in the neutral component of the state, i.e. c n = δ n1 , and 7 In the literature the Sommerfeld enhancement is sometimes quantified as S ≡ |ψ(0)| 2 /|ψ (0) (0)| 2 , the ratio of the potential-modified wavefunction at the origin to the wavefunction of a plane wave, ψ (0) = e ikz .
The presence of the mass shift, δM , causes the charged component of the state to decay exponentially at large r. A generic wavefunction in a spherical potential can be decomposed in terms of a radial wavefunction R kl (r) and Legendre polynomials as, We wish to determine the unknown coefficient, A n , by matching to Eq. 37 as r → ∞. For this, we expand the incoming plane wave into partial waves and also use the general form of the radial wavefunction for central potential scattering at long distances, where b n is a constant and δ n is the "phase shift." Strictly speaking, only the radial component of ψ 1 asymptotes to this form. As mentioned above, the radial component of ψ 2 decays exponentially. We can analytically continue k 2 and δ 2 though, to complex values to keep the decaying solution in the form of Eq. 40. Formally, δ 2 diverges to cancel the exponentially growing mode for complex k 2 , but as we will see, we never need to input its value to determine Sommerfeld enhancement.
To proceed with matching, it is useful to include the other, regular, linearly-independent solution to the Schrödinger equation (generically an N -component state has N regular and N irregular solutions, that latter will be useful for our numerical analysis as described below).
Since we were already considering the case of an incoming neutral state in Eq. 37, c n = δ n1 , we can take the orthogonal solution to be that of an incoming charged state, c n = δ n2 . Thus, anything with a n index becomes a tensor with indices, n, j (except k n , since it is kinematic and thus invariant across solutions), where j determines whether the incoming plane wave is neutral (j = 1) or charged. For c, we therefore get c nj = δ nj , and j iterates between an incoming χ 0 χ 0 state, what we were already considering, or an incoming χ + χ − . We can now define matrices, C nj ≡ c nj /k n and (M ) nj ≡ e −i(δ ) nj (b ) nj . This allows us to find the matching coefficient in Eq. 39, where the first column of A gives the coefficients for our state of interest.
In the opposite limit, R kl (r → 0) ∼ r , and since our only concern is the wavefunction at the origin, we need only keep track of the S-wave component. As a further simplification, we define χ k ≡ rR k, =0 , such that the Schrödinger equation reduces to Combined with our matching coefficient, this gives a matrix of Sommerfeld enhancement factors, where j labels the external state and n determines whether the 00 or ± component undergoes perturbative annihilation.
As setup currently though, we are on the hook for calculating the phase shifts, δ nj , and amplitudes, b nj , of the asymptotic solutions. To remove this difficulty and any subtleties about the decaying nature of the charged-component term, it is useful build the Wronskian with the other two linearly-independent solutions of Eq. 42, those that give irregular radial wavefunctions. Denoting them as (χ ) nj (r), we demand Dropping , since we are only interested in the S-wave case, W =χ χ −χ χ. W is invariant with respect to r, and so we evaluate it at 0 and ∞, where we have used that the C matrix is diagonal, with elements k −1 m . Finally, we can get a simpler expression for Sommerfeld enhancement, since equating W (0) = W (∞) gives However, the LHS is exactly what we obtained in Eq. 43. Thus, for numerical analysis, we only need to find the neutral component of the wavefunction at infinity,χ 11 (r → ∞) and Furthermore, sinceχ 1i (r → ∞) ∝ e ikr we can extract the phase part, defining ξ i ≡χ 1i e −ikr such that ξ (∞) = 0, giving an even simpler numerical implementation.
We take the reduced Schrödinger equation, Eq. 42, and solve for ξ i (≡χ 1i e −ikr ) andχ 2i with the potential in Eq. 35 and boundary conditions at the origin given in 47. Defining x ≡ kr, where k = M χ v/2, this gives the following coupled equations, We solve them numerically, obtaining quantitative agreement with the earlier literature [7,8]. 8 . There are nonetheless a couple provisos to the analysis. Firstly, Eq. 48 requires us to input a WIMP velocity. In accord with previous references, we have chosen v = 10 −3 .
Scanning over a range of velocities, we found insensitivity to the precise number as long as it was below M W /M χ and we were not in the immediate vicinity of the one of the resonance peaks in Fig. 6. However, since the annihilation rate is so large at these peaks, a WIMP mass at these values is sufficiently ruled out to make this detail beyond our scope. Additionally, the formal boundary conditions demandχ 2i (X) = 0 and calculate the Sommerfeld enhancement from ξ i (X), with X → ∞, but in practice we must take finite X. For our M χ range of interest, the longest decay length for the charged state is determined by the δM plus Coulomb part of the potential, and is given by x dec. = 1/ 2M χ δM/k 2 − 1. We find that taking X = O(10)x dec. leads to numerical stability.
Looking at the results of our Sommerfeld analysis in Fig. 6, we see the expected resonance structure, with values of |ψ(0)| 2 that exceed 10 4 . The Sommerfeld factors |ψ 2 (0)| 2 and |ψ 1 (0)| 2 are comparable throughout our range and maintain the modest hierarchy of the perturbative charged-state annihilation due to f + exceeding f − and its contribution having a larger numerical prefactor in Eq. 28. former is promotional to |ψ ± (0)| 2 and the latter to |ψ 00 (0)| 2 , as shown in Eq. 36.

VI. DARK MATTER CONSTRAINTS AND CONCLUSION
Having calculated tree level matching, LL resummation, and computed the Sommerfeld enhancement numerically, we can now evaluate the differential cross section for χ 0 χ 0 → γ + X, given in Eq. 28. We plot this in Fig. 7, where we have digitized the HESS limits given [24]. We note that in contrast to those groups that performed an exclusive two-body FIG. 7. Annihilation cross section to γ + X. Exclusion taken from [24], assuming an NFW profile.
calculation, [23,24], we find the effect of higher order correction to be very modest.
where for each value we have included Sommerfeld enhancement, and "NLO-fixed" includes only those one-loop effects that are resummed by our LL operator running. Thus, at this value for M χ , the leading corrections shift the semi-inclusive annihilation by just a few percent. Comparing directly to [24], which also investigated annihilation of a triplet fermion, they find that at 3 TeV, higher order effects lead to a ∼50% reduction in the exclusive rate.
This difference is to be expected given the distinct difference in our choice of observables.
From Eq. 28 and Fig. 6, we see that the leading contribution by a factor of few to our rate in our range of interest comes from perturbative χ + χ − annihilation and is proportional to |ψ ± | 2 ] . Thus f + , which has a tree level value of 2, drops to 1 in the limit M χ , E γ → ∞, meaning the |ψ ± | 2 term decreases by 50% only in this infinite limit. Furthermore, in this asymptotic case, f − → 1 from its tree-level value of 0, so the |ψ 00 | 2 and ψ ± ψ 00 terms can boost overall rate above this 50% reduction. Comparing to the LL results in Eqs. 13-16 of [24], their Sudakov factor goes like exp(− α W π log 2 ( M W Eγ )). Thus, in the limit of infinite DM mass, their rate drops to 0. This is expected from the general result that exclusive rates vanish in the limit of infinite energy.
The limit from HESS in Fig. 7 shows that the thermal relic wino, M χ ≈ 3 TeV is ruled out by more than an order of magnitude. Unfortunately, the astrophysical uncertainties in the halo profile are sufficient to evade an excess of even this size. This is because the flux of photons measured by the HESS experiment is proportional to the "J-factor", where ρ loc is the local density, R = 8.5 kpc is the distance to the galactic center, and s is the line of sight distance to the experiment, where r = s 2 + R 2 − 2sR cos θ. Discussions on the ability of different halo models to evade constraints can be found in the earlier papers that found the wino to be in tension with HESS [4,7]. The exclusion curve we have taken from [24] assumes an NFW profile [32] with a local density, ρ loc = 0.4 GeV/cm 3 [33], and r s = 20 kpc [34], This is a cusped profile, diverging as 1/r toward the galactic center. 9 In the discussion that follows, we fix the local density ρ loc = 0.4 GeV/cm 3 , but we will change the functional form of the distribution along with a possible core radius. It is possible though, that the local density could lie somewhere in the range of 0.2-0.6 GeV/cm 3 [34]. The cusped vs. cored (where the distribution flattens out at some distance) debate on the nature of the DM halo is an old one that experiment is far from resolving. A recent observational analysis found good fits both for NFW profiles similar to the one in our constraint and for ones with relatively large, ∼10 kpc cores [36]. Looking at simulation, DM-only models generically yield cusped distributions [37]. However, observations of dwarf spheroidal galaxies found evidence for cored profiles, which was subsequently found in numerical models that included effects from baryons [38]. It is thought that supernovas near the galactic center may eject enough mass to flatten the DM distributions. The question is whether sufficient mass could be ejected in a much more massive, (100 − 1000)× larger, Milky Way-like galaxy, or if the larger baryon density near the galactic center results in an even more cusped distribution.
Simulations of Milky Way-like galaxies that include baryons continue to show a preference for cusped distributions, at least into distances ∼1 kpc [38][39][40]. One can ask, how much coring is needed to save the wino, given our LL-resummed annihilation rate? For an NFW profile that becomes constant below a certain radius, (r/rs)(1+r/rs) 2 r > r c ρ 0 (rc/rs)(1+rc/rs) 2 r ≤ r c in Fig. 9 we plot the value of the core radius, r c , needed to make our semi-inclusive rate calculation consistent with the limit from HESS. For the thermal relic wino mass, M χ = 3 TeV, if we consider the Burkert profile [41], where r c again gives the core radius, and a cutoff-NFW distribution, we find that we need cores of 4-4.5 kpc and 1-1.5 kpc, respectively, to avoid the bounds. Both distributions are 9 Typically, there will be deviations from strict spherical symmetry in the halo (axial, triaxial), and these can effect J at the 10-20% level [35]. This uncertainty though, is far below the orders of magnitude shift in J one can obtain by changing the profile shape between cusped and a large core.
well within the region allowed by observation [36], and the wino in a cutoff-NFW profile (but not Burkert) is consistent with simulation, as well.
The observation of wino dark matter near the thermal relic mass of 3 TeV would point to the existence of a nontrivial amount of coring in the halo of the galaxy which would require an explanation. Of course, there are other possible ways to evade the HESS constraints, even if the profile were nearly NFW. There is the possibility that the lightest neutralino may not be a pure wino. For example, a thermal relic higgsino is far from constrained, and thus admixtures between these states could certainly be allowed [4]. Sticking with the pure wino, if there were some non-thermal mechanism for its production, then the limit at values other than 3 TeV would be relevant, and M χ could be in one of the allowed regions shown in Fig. 7. Alternatively, whether or not its production were thermal, the wino could make up just a fraction of the dark matter, and thus much of parameter space would remain open, as shown in Fig. 8. With the theoretical uncertainty on its annihilation rate now under control at the O(1%) level, 10 the discovery of a wino at future indirect detection experiments, such as CTA [42], could give us important windows into further open questions such as the halo 10 It would be an interesting exercise to extend this analysis to NLL. We have computed the running of our Wilson coefficients from the one-loop cusp anomalous dimensions. One would also need one-loop noncusp, two-loop cusp, and the β-function running of α W . These were included in the exclusive-observable calculations of [23,24]. Additionally, the one-loop running of our fragmentation functions, Eq. 12, is needed.
FIG. 9. The amount of coring required for the wino to become viable with respect to the HESS constraint shown in Fig. 7 for the cutoff-NFW profile (Eq. 52). The three curves display the effect of variation in the local dark matter density.
distribution, cosmological history of DM production, and the presence of multi-component dark matter.
This is an approximation to a full, quantum field theoretic treatment, but correctly sums all ladder-exchange diagrams [8,27,28]. To obtain Eq. A1 from the relativistic, electroweak lagrangian for an SU(2) triplet, we must define the two-body states evolved by the hamiltonian containing V (r) in terms of the field theory creation and annihilation operators. After finding the four-fermion interaction in the nonrelativistic theory by matching to potential boson exchange, we can then find V (r) by calculating how the interaction modifies our properly normalized, two-body states. This allows us to write it in a form which makes no reference to fields, but evolves the state directly as a function of r.
We start with the standard lagrangian for a nonrelativistic field theory, which can be obtained for fermions by projecting to the single-particle sector of the Hilbert space (therefore decoupling the particle and antiparticle fields), pulling out a factor of e −iMχt and then integrating out the small components of the fermion 4-spinor, leaving a field theory of 2spinors, We have adopted the convention that χ 0 and χ − in Eq. A2 destroy their respective particles, while χ + is a creation field. Explicitly, we have where p 0 = p 2 /2M χ in the Fourier factor and η s is a Pauli spinor with η 1 = (1 0), η 2 = (0 1), and η −1 = η 2 , η −2 = η 1 . Below, we use arrows to denote the spin. Additionally, since χ 0 is a Majorana fermion, it is useful to define the following field, which also annihilates χ 0 particle, Matching to the full theory we consider the exchange of a potential mode with propagator i/(∆p 2 + M 2 ), where ∆p = p − p is the difference in spatial momentum of the incoming and outgoing fermions and M is the mediator mass (Fig. 10). In terms of our nonrelativistic fields, and projecting to spin-singlets for our fermion bilinears via a Fierz transformation, we get which agrees with [8]. To pass from nonrelativistic field theory to quantum mechanics, we need to construct our two-particle states. Since the subtlety of obtaining Eq. A1 concerns the overall factors, we must pay special attention to the overall normalization. For a state, |B( p) s , of mass M B and spin s, we have that B( p ) r |B( p) s = 2M B (2π) 3 δ (3) ( p − p)δ rs .
Since we work in the CM frame of our two particles, there is no total momentum, and we therefore get the normalization factor 2M B (2π) 3 δ (3) (0)δ rs . Starting with our charged state, we satisfy Eq. A6 if we normalized it as where ψ(q) is the relative momentum wavefunction, whose Fourier transform we will evolve with V (r). The calculation for the neutral, two-particle state is similar, but since χ 0 is its own antiparticle, we get twice as many contractions, and the overall factor is therefore √ 2 smaller, Finally, extracting the field theory potential from S pot'l in Eq. A5, we can check how the interacting hamiltonian acts on states |N , |C . As an example, we take the off-diagonal interaction that converts a neutral state to a charged state and count the factors of 2. The potential in S pot'l contains an overall 1/2, and there are 4 possible contractions of the d operators in the interaction with the d † operators in the state. Finally the state |N has an overall normalization of 2M χ and a 1/ √ 2 from the spin-singlet normalization. Thus, we get 2× the creation operator structure of the |C state. However, its normalization is an overall 4M χ times the 1/ √ 2 from being in the spin singlet, for an overall √ 2. Thus, the potential that evolves |N into |C requires an overall √ 2, as we see in the "12" entry of Eq. (A1). Factors for the other entries follow similarly. We see that the nontrivial √ 2 in the off-diagonal entry originates from the mismatch in normalization between the neutralinos, where the two-body state contains identical particles, and the chargino, where the two particles are distinct.