Sommerfeld effect for continuum gamma-ray spectra from Dark Matter annihilation

We present a calculation of the continuum part of the gamma-ray spectra resulting from Dark Matter annihilation in the framework of the MSSM taking into account Sommerfeld effects. Concentrating on pure wino and pure higgsino scenarios we compare our calculation to existing work and explore the numerical impact of the features not captured by previous approximative descriptions. We find that, in particular for large neutralino masses, when the Sommerfeld enhancement is very large, chargino-antichargino annihilation processes, which have not been considered before, lead to sizable differences with respect to existing calculations. In scenarios with neutralinos in the intermediate-mass range, we find that the role of the charginos is crucial in the endpoint regime. Our calculation provides the currently most accurate prediction for the continuum gamma-ray spectra.


Introduction
Understanding the nature and origin of Dark Matter (DM) represents one of the great challenges of contemporary research in cosmology, astro-and particle physics.While there is little doubt in the community on the very existence of such a form of matter that does not interact electromagnetically with ordinary matter, the Standard Model of particle physics (SM) does not offer an obvious candidate particle to account for it.
Promising DM candidates are provided by models of physics beyond the SM (BSM) featuring weakly interacting massive particles (WIMPs) with masses in the range of a few hundred GeV to a few TeV.It is assumed that after a period of equilibrium between thermal WIMP production and annihilation in the early Universe a freeze-out of their number density occurred when the expansion rate of the Universe became larger than the relevant interaction rate.The measured abundance of DM is compatible with the existence of a WIMP at the electroweak (EW) mass scale.
A natural framework for WIMPs is provided by supersymmetric theories that complement the particle spectrum of the SM by so-called superpartners with spin differing by one half.The breaking of supersymmetry (SUSY) allows these particles to acquire masses different from their SM partners.In many SUSY models the lightest stable particle (LSP) represents a viable DM candidate.A prime example of such a model is constituted by the Minimal Supersymmetric Standard Model (MSSM).The MSSM contains a minimum number of SUSY particles, i.e. one superpartner for each SM particle apart from the Higgs boson.The Higgs sector has to be extended, resulting in five physical bosonic Higgs states and the same number of fermionic higgsinos.Form a theoretical point of view, the MSSM is appealing not only because of its underlying symmetry structure, but also by ensuring the unification of gauge couplings, naturalness, and providing a DM candidate.In most realizations of the MSSM the lightest neutralino, which effectively is a mixture of bino, wino and higgsino components, is considered stable and thus constitutes a viable WIMP DM candidate.
Great hopes for the discovery of SUSY particles were pinned on the CERN Large Hadron Collider (LHC) with the highest collision energies ever achieved on Earth by a proton accelerator.Dedicated searches by the ATLAS and CMS experiments put limits on the MSSM parameter space, but still leave room for the existence of SUSY particles, in particular in socalled split scenarios where the LSP is much heavier than the scale of EW symmetry breaking [1,2,3].In such scenarios neutralino mixing effects are typically suppressed and the DM candidate is part of a "pure" EW multiplet [4].Such scenarios are appropriately referred to as higgsino or wino type.Due to the large neutralino masses these DM models predict they remain elusive to searches at colliders and direct detection experiments.
Depending on the underlying model, DM candidates are supposed to be produced in pairs or in association with accompanying particles in collisions of SM particles [5,6].If the masses of the DM particles are large and/or their interactions with SM particles feeble, the corresponding production rates are low.The same limitation applies to direct detection experiments that aim to identify the recoil of a DM particle off a nuclear target [7,8,9,10].
Promising alternatives for DM searches in such scenarios are constituted by indirect detection strategies aiming at identifying annihilation signatures of DM particles [11].In particular, when such annihilation processes are accompanied by photon emission the resulting gamma ray spectra feature a very characteristic peak structure at the kinematical endpoint.The shape and intensity of such spectra can be heavily influenced by the so-called Sommerfeld effect.This, in turn, can be exploited to overcome severe limitations of indirect searches due to large uncertainties in the distribution of DM in the inner galaxy and omni-present astrophysical backgrounds.The Sommerfeld effect is ubiquitous in annihilation processes involving non-relativistic particles that can exhibit long-range interactions.This phenomenon applies to MSSM neutralinos, where the long-range interactions are mediated by electroweak bosons.
A plethora of experiments and astrophysical observations has been devised to make use of this detection strategy.Particularly interesting for indirect neutralino searches are various observatories.These include the space telescope Fermi-LAT [12] which is suitable for neutralino searches in the mass range of O(1) − O(100) GeV.Additional information is coming from the currently operating Imaging Air Cherenkov Telescopes H.E.S.S. [13], VERITAS [14], MAGIC [15], along with their next-generation counterparts CTA [16] and LHAASO [17].The water Cherenkov telescope HAWC [18] further enriches this list.These Cherenkov telescopes are particularly suited for neutralino searches with somewhat heavier masses in the range of O(0.1) − O(100) TeV.
From the theoretical point of view, exploiting the full potential of this search strategy requires a quantitative understanding of the Sommerfeld effect.In the context of the full MSSM, work to that effect has been limited.In particular, Refs.[19,22] are, to our knowledge, the only papers including the Sommerfeld enhancement in their MSSM analyses.In these studies, however, the incorporation of the Sommerfeld effect is achieved through the matching of exclusive 2-to-2 Sommerfeld-resummed neutralino-annihilation cross sections.This approach has some caveats and, as we will argue, cannot account for several pivotal phenomenological aspects that can become crucial for future analyses.
We aim to fill that gap with this work.In particular, we compute all neutralino and chargino annihilation cross sections into three-body final states that are relevant for obtaining the continuum gamma-ray spectra.We obtain these annihilation cross sections in analytical form for arbitrary spin and helicity combinations of the final-state particles.
While our results are valid for generic MSSM parameter sets, for our numerical discussion we focus on the pure wino and higgsino limits.We find that the impact of the Sommerfeld effect on the continuum spectrum is sizable even when the neutralino mass is of the order of a few hundred GeV.Our results show that, besides the Sommerfeld effect being very large [32], the chargino contribution dominates for the intermediate to the very high energy part of the gamma-ray spectrum.This aspect has not been captured by previous computations.
The paper is structured as follows: In Sec. 2 we briefly review the basic theoretical aspects of indirect DM detection using gamma rays.We then move on to the discussion of the Sommerfeld effect in the context of SUSY and the methods we used for the calculation of annihilation cross sections in Sec. 3. In Sec. 4 we discuss our numerical results in the context of pure wino and higgsino scenarios, and we then conclude in Sec. 5. Conventions and some technical aspects of our work are discussed in App.A and App.B.
where J obs is the astrophysical "J" factor [33] for a given observed region and ⟨d(σv)/dE γ ⟩ is the velocity-averaged annihilation cross section of two DM candidate particles χ of mass m χ into gamma rays.With some rare exceptions the J factors are independent of the gamma-ray energy E γ , and are thus irrelevant for the description of the spectral properties of the DM gamma-ray signals, which are the focus of this work.We refer the reader interested in specific J factors to Refs.[34,35,36] for the Milky Way and to Refs.[37,38,39] for other astrophysical environments.
Here, however, we concentrate on the DM annihilation cross section d(σv)/dE γ with the DM candidate particle being constituted by the lightest neutralino χ 0 1 of the MSSM.This annihilation cross section has a kinematic endpoint at E γ ≃ m χ1 .At this energy, the spectrum features a quasi-monochromatic gamma-ray line with a (natural) broadening of O(v 2 ), where v is the average speed of the neutralinos in the DM halos of interest (v ≈ 10 −3 in the Milky Way).The monochromatic nature of the line is due to the two-body kinematics of the neutralinoannihilation process into photons, χ 0 1 χ 0 1 → γγ.On top of the gamma-ray spectral line, the endpoint spectrum of neutralino annihilation involves a Z resonance (from the χ 0 1 χ 0 1 → γZ * process) with a natural width of Γ Z /m Z ∼ 0.03.Both the γγ and the γZ contributions at the endpoint of the spectrum are, in principle, loop suppressed.However, the narrow width of the Z resonance and the significant influence of long-range interactions between the neutralinos and the charginos due to the Sommerfeld effect make these features highly intriguing in terms of detectability.
The remaining part of the spectrum, the continuum, is generated by the sum of all processes where neutralinos annihilate into a gamma-ray photon in association with a multiparticle con- Illustrative Feynman diagram for the continuum contribution to the gammaray emission spectrum produced by neutralino annihilation.Red lines depict final-state photon radiation, blue lines soft/collinear electroweak radiation, and purple lines internal bremsstrahlung.figuration "X", i. e. χ 0 1 χ 0 1 → γ + X.At leading order (LO) in the electroweak coupling, X consists of two SM particles.For clarity, we use a superscript in order to differentiate such a two-body state X (2) from the more complex sub-states X that can occur in general, and that are illustrated by Fig. 1.
In many simulations, such continuum configurations are approximated by matching fixedorder computations for two-particle production processes to parton-shower programs such as Pythia [40] or Herwig [41].In this approach gamma-ray spectra from WIMP annihilation are obtained using the parton-shower approximation formula where, given a specific WIMP model, the coefficients (σv) X (2) are the (tree-level) cross sections for the annihilation of two DM candidate particles into two SM particles (e.g.X (2) = b b, τ + τ − , W + W − , . ..), while the model-independent functions dN MC X (2) →γ /dE γ for the gammaray emission off the SM particles are obtained from the aforementioned Monte-Carlo event generators. 2 These functions are available in specialized software packages such as DarkSUSY [42], PPPC [43] or MicrOMEGAs [44].
As a diagrammatic visualization, in Fig. 1 we show an example of a Feynman diagram for the continuum contribution to the gamma-ray emission spectrum produced by neutralino annihilation.From that particular diagram only those photon lines that are emitted from final-state legs are accounted for by the early implementations of Eq. ( 2) (e. g.DarkSUSY versions before v5 [45]).Such contributions are termed final-state radiation (FSR).Soft/collinear electroweak radiation effects have also been included in the literature (see, e.g., Ref. [46]).In particular, the fragmentation functions provided by the PPPC code include these corrections.Newer versions of DarkSUSY instead provide a more complete picture for the neutralino-annihilation photon spectrum at the expense of losing the model-independence of Eq. (2).In addition to pure FSR calculations they capture potentially dominant processes such as the so-called internal bremsstrahlung (IB) [47,48,49,50,51], which is absent in the PPPC approach [43].This is achieved by computing the full fixed-order χ 0 1 χ 0 1 → γ+X (2) annihilation cross section for a given state X (2) (X (2) = W + W − in the sample diagram of Fig. 1) and matching it to a parton-shower program while carefully subtracting redundant terms in order to avoid double counting, see e. g.Ref. [50].In particular, IB becomes important in those cases where the otherwise helicitysuppressed annihilation of non-relativistic Majorana particles into a particle-antiparticle pair of light fermions becomes sizable once radiation effects are accounted for.It has been even observed that in some WIMP models, IB gives the dominant contribution to the gamma-ray spectrum in the medium-to-high energy regime (see Refs. [52,53,54] for more details).

Sommerfeld effect
The previous discussion evidences that great progress have been achieved in understanding the continuum part of the gamma-ray spectrum from neutralino annihilation.By combining fixedorder computations with Monte-Carlo event generators (see Eq. ( 2)), a relatively adequate theoretical picture of the WIMP gamma-ray spectrum can be obtained.This picture, however, fails at capturing virtual (loop) effects which, as we will see below, can become crucial.In particular, multi-loop corrections such as the one depicted in Fig. 2 can have an enormous impact on the annihilation cross sections of heavy neutralinos [55,32,56], resulting in enhancement factors of several orders of magnitude.
The primary reasons for this phenomenon are (A) the non-relativistic nature of the initial state consisting of two neutralinos, (B) t-channel interactions between them mediated by the gauge and Higgs bosons of the theory (see sketch in Fig. 2), and (C) the fact that these exchange bosons are lighter than the annihilating neutralinos.These features make a quantummechanical approach in terms of static potentials the most appropriate in order to account for the self-interactions of the neutralinos prior to their annihilation.In more formal terms, the computation of Feynman diagrams such as the one depicted in Fig. 2 will inevitably yield terms that are (parametrically) 3 , where n is the number of loops in the given ladder-like diagram and α denotes the fine structure constant.Therefore, since we assume that m W ≪ m χ , the contribution of each of these diagrams is large and cannot be neglected.Rather all such terms need to be systematically resummed.
The resummation of these diagrams can be performed consistently in the context of nonrelativistic effective field theories (NREFTs) [57] (in the context of MSSM neutralino annihilations, see e. g.Refs.[58,59,60]) independently of the exclusiveness of the final state.For the  3) for the representative particular case of neutralino annihilation into photons, the resulting annihilation cross section is given by [61,62] where S IJ is the matrix of the so-called Sommerfeld factors, and we will refer to d( σv)/ dE γ as the annihilation matrix for the Sommerfeld-corrected χ 0 1 χ 0 1 → γ + X process.The pair indices I, J denote all possible neutral combinations of neutralino-neutralino and chargino-antichargino pairs within the MSSM [60].There are 14 such combinations: We compactly denote each pair index K by nested particle indices distinguishing between neutralino and chargino states using either round or angular brackets, i.e.K = {( îj)}, where ( îj) = ( 11), ( 12), . . .for the four neutralinos (i, j = 1, . . ., 4) and K = {⟨xȳ⟩}, with ⟨xȳ⟩ = ⟨1 1⟩, ⟨1 2⟩, . . .for the two charginos (x, y = 1, 2).Note that while the neutralino states satisfy ( ĵi) = ( îj), the chargino states are not exchange symmetric (⟨yx⟩ ̸ = ⟨xȳ⟩).Figure 3 captures the essence of Eq. (3) in the representative χ 0 1 χ 0 1 → γ + f f process.In that particular example, the corresponding (interference) term of the equation with I = ⟨xȳ⟩ and J = ( îj) is shown.The diagram explicitly illustrates that both the Sommerfeld factors and the annihilation matrix elements can be expressed as products of quantities corresponding to the intermediate states I and J, respectively.In the calculation of the annihilation matrix elements, the amplitudes of the χ + x χ − y → γ + f f and the (complex conjugated) χ 0 i χ 0 j → γ + f f processes enter.For the computation of the S IJ factors quantum mechanical wave functions evaluated at the spatial origin have to be multiplied with their complex conjugates.These wavefunctions describe the non relativistic χ 0 i χ 0 j transitions as we briefly discuss in the following.

Computation of the Sommerfeld factors
Up to corrections of O(v 2 ) and O(m 2 W /m 2 χ ), the resummation of ladder-like diagrams such as the ones shown in Fig. 2 requires solving the following matrix Schrödinger equation 4 for the multicomponent function u IK (r) [60] The Sommerfeld factors S IJ in Eq. ( 3) are then given by where ) and M I is the sum of the masses of the non-relativistic "I" state, e. g.
The potential matrix in terms of the "coupling matrices" αB IJ is given by while the "mass-splitting" matrix ∆M in the basis of Eq. ( 4) reads where α is the fine-structure constant.The remaining coupling matrices depend on the MSSM parameters through the neutralino and chargino mixing matrices defined in Eqs. ( 24)-(25) of App. A. These had been already obtained in Ref. [60], and we re-derived them using the conventions of Ref. [63] that we use throughout this work.In the Feynman gauge the entries of these coupling matrices read for the vector mediators Z and W ± and for the scalar and pseudoscalar mediators, where S collectively labels all physical neutral Higgs bosons (S = h, H, A 0 ) in the MSSM, and the would-be Goldstone bosons G Z , G W associated with the Z and W bosons.
Explicit expressions for all the coefficients v B IJ , a B IJ and s B IJ can be found in App. A. We obtained these by using the Mathematica [64] packages FeynArts [65] and FormCalc [66,67] with the MSSM model file [63].At the relevant perturbative order, our potential agrees with the corresponding results of Ref. [60].

Computation of the annihilation matrices
In the previous section we briefly reviewed the NREFT methods that are necessary in order to properly incorporate the Sommerfeld effect in the continuum spectrum prediction.In particular, Eq. ( 3) provides us with an elegant prescription on how to deal with this problem.Given the Sommerfeld coefficients S IJ which can be computed by solving a system of Schrödinger equations, we need to determine the [d( σv)/dE γ ] IJ functions for every possible combination I, J.The IJ element of the annihilation matrix is defined by where the identical-particle index id(K) = δ ij for neutralino pairs [K = ( îj)] or 0 for charginoantichargino pairs [K = ⟨xȳ⟩].dΠ ′ γ+X is the phase-space integration element for the γ +X final state, and A (l,s)=(0,0) K→γ+X is the amplitude for the s-wave annihilation5 of the two-particle state K into a γ + X state.It is implicitly understood that the sum over X includes all possible spin and helicity combinations of the particles that are produced in association with the gamma ray.
In this article we compute the annihilation matrices of Eq. ( 17) for the production of a photon in association with any possible combination of two SM particles X = X (2) .Concretely, in the MSSM the only possible combinations are6 where S and S ′ are shorthand for two different neutral Higgs scalars in the MSSM and f f denotes any fermion anti-fermion pair within the SM.There are 14 × (14 + 1)/2 = 105 independent symmetric combinations of the initial-state indices (see Eq. ( 4)).Thus, for each one of the aforementioned combinations of X (2) 105 independent matrix elements have to be computed.Out of these matrix elements, in the literature only one is available [50], corresponding to the neutralino-neutralino annihilation cross section with I = J = ( 11).In this work we compute the annihilation matrices for the remaining 104 combinations of I, J.
State-of-the-art software packages such as FeynArts [65] or FormCalc [66,67] are capable of computing differential cross sections in analytical form for annihilation processes within complicated models such as the MSSM.However, the problem at hand requires some extra processing for the s-wave projection of the I, J states and for computing the interference of amplitudes with different initial states.
We addressed these issues by obtaining raw amplitudes for the 2-to-3 scattering processes χ 0 i χ 0 j → γ +X (2) of all possible X (2) states of Eq. ( 18) with arbitrary spin and helicity combinations using FeynArts 3.11/FormCalc 9.8 and the built-in MSSM model file of Ref. [63].We then processed these amplitudes in the framework of Mathematica [64] to obtain the desired interference contributions.

Annihilation matrices in the FSR approximation
The calculation of the annihilation matrix d( σv)/ dE γ presented above is exact at LO in the electroweak couplings.Thus, when multiplied with the Sommerfeld factors and matched with parton-shower simulations, the resulting prediction offers the most accurate description of the continuum gamma-ray spectrum from annihilating neutralinos up to the present day.However, an approximate treatment that exploits the fact that the computation of the Sommerfeld factors is independent of the exclusiveness of the final state in the annihilation process, is also possible and has already been considered in Refs.[19,22].In that approach, Eq. ( 2) generalizes to dN MC where, in analogy to Eq. ( 3), the Sommerfeld-corrected 2-to-2 annihilation cross section is given by (σv and, as just argued, the fragmentation functions remain the same as in Eq. ( 2).Eq. ( 19) is sensitive to the different polarizations of the 2-body final states.To emphasize this feature, the polarization states can be indicated explicitly as where the subscripts "T " and "⊙" refer to their transverse and longitudinal polarization components of the gauge bosons and the "L" and "R" subscripts denote the left-and right-handed chirality of the fermions, respectively.While the unpolarized cross sections for these 2-body states are already known [58], we obtain the polarized ones for the first time in this work.

Numerical results
In light of the inherent complexity of the MSSM, a thorough exploration of our calculations requires a dedicated study.In this section, we focus on two limiting scenarios of the MSSM: the pure wino and the pure higgsino case.These models have been investigated intensively in the last several years (see e. g. [19,55,68,69,70,71,72]).We note that even in the most generic MSSM parameter sets, the Sommerfeld effect has a significant impact on the associated indirect detection signals.
The most general form of the softly broken MSSM introduces 105 parameters in addition to those present in the SM [73].These are interdependent given a specific supersymmetry breaking scenario.Current experimental efforts to search for SUSY, however, concentrate on constrained versions of the MSSM such as the 18-parameter phenomenological MSSM (pMSSM) [74] or the 4-parameter constrained MSSM (CMSSM) (see, e.g., Refs.[75,23]).The pMSSM, for example, treats the mass terms for gauginos (M 1 , M 2 , M 3 ), higgsino (±|µ|), sfermions and trilinear couplings, as well as the ratio of the vacuum-expectation values of the two Higgs doublets (tan β), as independent parameters.
In the pure wino and higgsino limits, the sfermion masses as well as the trilinear couplings, the gluino mass parameter M 3 , and tan β are assumed to be infinitely large.In the pure wino limit, in addition to integrating out the sfermions it is assumed that M 1 → ∞ and |µ| → ∞, leaving M 2 as the only relevant MSSM parameter.In the spirit of the minimal DM models discussed in Ref. [68] (see also Refs.[76,77] for a more recent study on minimal DM) the wino can be visualized as an SU(2) Majorana triplet which, after electroweak symmetry breaking, gives rise to one chargino and one neutralino degree of freedom.Similarly, the pure higgsino limit can be constructed by introducing an SU(2) Dirac doublet that results in two neutralinos and one chargino.Note that the pure wino and higgsino scenarios do not provide satisfactory solutions to the hierarchy problem and the unification of gauge couplings, which were among the primary motivations for introducing SUSY.Moreover, for sub-TeV neutralino masses within the standard freeze-out hypothesis they predict more DM than observed [4].Nonetheless, these limiting scenarios serve as highly valuable toy models.For instance, the dimension of the Sommerfeld and the annihilation matrices in Eq. ( 3) is reduced to two for the pure wino and to three for the pure higgsino model, respectively, thus significantly reducing the complexity of calculations and analyses as compared to the full MSSM.
For the numerical results shown below as electroweak input parameters we choose the mass of the Z boson,m Z = 91.1876GeV, the Fermi constant, G µ = 1.1663788×10 −5 GeV −2 , and the fine structure constant, α 0 = 1/137.035999180as quoted by the Particle Data Group [78].From these input parameters, the mass of the W boson, m W = 80.360 GeV, and the Weinberg angle sin 2 θ W = 0.22338 are obtained using the electroweak relations of Ref. [78].In our numerical implementations we use the running QED coupling evaluated at the scale m Z , α = 1/128.93.The masses of the u, d, s quarks, the electron and muon, and of all neutrinos are neglected.For the c, b, t quarks, the tau lepton, and the Higgs boson we use the masses quoted in Ref. [78]: For both, wino-like and higgsino-like models, we use the one-loop expressions given in Refs.[68,79] for the mass splitting between the chargino and the lightest neutralino, m χ ± 1 − m χ 0 1 .For the mass splitting between the two neutralinos in the pure higgsino scenario we use m χ 0 2 − m χ 0 1 = 20 MeV in all our examples, as done in Ref. [80].The numerical results that are reported below involve comparisons between different prescriptions and approximations.For simplicity we introduce the following acronyms: • 'fixed noSE': Neutralino-pair annihilation into a γ +X (2) final state at tree-level with no Sommerfeld resummation at fixed order (O(α 3 )), as provided in Ref. [81].Generally, X (2)  denotes a two-particle state.In the pure wino and higgsino models, only • 'FSR SE': Sommerfeld-corrected neutralino-pair annihilation into a γ + X (2) final state where only photon emission due to FSR in the soft/collinear approximation is taken into account.

SE
Figure 4: Gamma-ray spectrum of neutralino annihilation in the higgsino-like (black) and the wino-like (gray) scenarios with m χ = 150 GeV within our full calculation (solid lines) and the unresummed LO computation of Ref. [81] (dashed lines).
For the pure wino and higgsino scenarios, the only possible two-body final state is The spectrum that is obtained using this approximation is computed by using Eq. ( 19) with the LO fragmentation function of Ref. [46] which in this case assumes the particularly simple form • 'full SE': Neutralino-pair annihilation into a γ + X (2) final state including continuum contributions and Sommerfeld resummation effects.In contrast to the two previous approaches, this computation receives contributions from all the two-particle combinations of Eq .(18).

Scenarios with m χ = 150 GeV
In Fig. 4 we show the gamma-ray spectra of neutralino annihilation in these two scenarios for our full calculation including the Sommerfeld effect and compare it to the unresummed LO calculation of Ref. [81].Note that the results of the 'FSR SE' approach are zero in the plot range due to the kinematic threshold in Eq. ( 21) at E γ ≈ m χ − m W ≈ 70 GeV.The comparison between the two remaining calculations illustrates impressively the relevance of the effects that have been neglected in the past, even when there are no large hierarchies between the masses of the lightest supersymmetric particles and the gauge bosons of the electroweak theory.The most striking peculiarity of Fig. 4 is the appearance of a resonance associated with a Z boson contribution which is omitted by the "naïve" fixed order computation.In particular, the fixed-order calculation has a clear cutoff from the kinematic threshold of the χ 0 1 χ 0 1 → W + W − γ process at the energy E γ = m χ − m 2 W /m χ (corresponding to about 107 GeV in the scenarios of Fig. 4).The Sommerfeld-corrected continuum spectra 'full SE' that we compute here account for this Z-boson contribution by the inclusion of Feynman diagrams like the one shown in Fig. 2. In such diagrams fermion pairs can be created resonantly from a Z-boson in association with a photon with an energy that is larger than the aforementioned cutoff.The resulting peak of the Z resonance is then at E γ = m χ − m 2 Z /(4m χ ), amounting to about 136 GeV in the scenarios of Fig. 4.
In addition to the Z resonance, the spectral region we uncovered in this work exhibits further interesting features.The almost imperceptible kinks around the E γ = m χ endpoint (150 GeV in the scenarios of this section) are due to the kinematic thresholds in the χ 0 1 χ 0 1 → f f γ processes, where f denotes any charged fermion of the SM with non-vanishing mass.For instance, the We note that, as expected, the impact of the Sommerfeld effect in the spectral region covered by the original fixed-order computation of Ref. [81], i.e. below the W + W − threshold, is marginal.As the m χ /m W hierarchies become larger, the numerical impact becomes more pronounced, as we shall discuss below.Independent of the neutralino mass, however, annihilation cross sections and Sommerfeld factors from wino-like neutralinos are generically larger than in the higgsino-like case, because the electroweak charges in the pure wino limit are a factor of two larger than those for the pure higgsino limit.For instance, the Born-level annihilation cross section of wino-like neutralinos into W + W − pairs is sixteen times larger than the corresponding cross section for higgsino-like neutralinos.

Scenarios with heavier neutralinos
Accounting for the Sommerfeld enhancement becomes more and more important as the mass of the neutralinos is increased.In Fig. 5 we demonstrate this numerically for higgsino-and wino-like scenarios with a lightest neutralino of mass 600 GeV and 2.4 TeV, respectively.
In each of these heavy neutralino scenarios, the predictions obtained within our full calculation show strong enhancements compared to the fixed-order computations.In the soft-photon part of the spectrum shown in Fig. 5 we observe how well the FSR predictions 'FSR SE' match our improved predictions 'full SE' and how this agreement improves, as expected, when m χ is increased.In particular, due to the model-independence of the splitting function of Eq. ( 21), the ratios of the Sommerfeld-resummed calculations 'FSR SE' to the fixed-order results 'fixed noSE' are roughly independent of E γ as long as E γ ≪ m χ .In the higgsino-like scenario these ratios amount to an enhancement by a factor of about 1.5 for m χ = 600 GeV and of about 4.6 for m χ = 2.4 TeV.In the wino-like scenario we find enhancement factors of ∼ 1.9 for m χ = 600 GeV and of 990 for m χ = 2.4 TeV.
This huge enhancement of the gamma-ray spectrum in the large-m χ pure wino scenario is due to the resonant nature of the Sommerfeld effect (see e. g.Fig. 1 of Ref. [32]).For particular neutralino masses the binding energy of a chargino-antichargino bound state is exactly zero, which ultimately results in a strong enhancement of the DM annihilation cross section.When varying the neutralino mass within our pure-wino scenario the enhancement factor reaches a  1: Elements of the Sommerfeld matrix in the pure-wino limit of the MSSM for the three neutralino masses considered in this section.
first resonance at m χ = 2.29 TeV. 7It then decreases and increases again, until the second resonance is reached at m χ = 8.83 TeV.Note that, as discussed e. g. in Ref. [83], depending on the SUSY parameters, the Sommerfeld resummation can lead to suppressed rather than enhanced annihilation rates.
As E γ increases, the situation becomes even more intriguing, requiring us to consider the relative impact of the various entries in Eq. ( 3) for unraveling the underlying dynamics.More specifically, in the pure wino limit Eq. ( 3) can be decomposed into three terms: where a numerical evaluation of the Sommerfeld factors yields the results listed in Tab. 1.In order to interpret the entries of the table, it is important to bear in mind that in the absence of the Sommerfeld effect S IJ = δ I( 11) δ J ( 11) .In the "low mass" scenario (m χ = 150 GeV), this condition is satisfied at the level of about 10 %.This, for instance, explains why in Fig. 4, the 'fixed noSE' and the 'full SE' computations are similar for E γ < 107 GeV.However, as noted before, the 'fixed noSE' calculation misses the Z resonance effect and all additional f f γ contributions resulting from the last term of Eq. (22).As the neutralino mass increases, the numerical significance of this term as well as of the interference term (second term in Eq. ( 22)) grows substantially.As shown below, the role of the charginos becomes more and more crucial as the neutralino mass is increased.Our 'full SE' calculation exhibits an interesting phenomenon at gamma-ray energies that are somewhat smaller than the threshold energy of the W + W − γ final state, E γ = m χ −m 2 W /m χ .Specifically, unlike the relatively "soft cutoff" behavior that is observed in the 'fixed noSE' computation, the 'full SE' prediction exhibits a subtle additional enhancement followed by a sharper decline near the threshold.This additional enhancement becomes more pronounced as the neutralino mass is increased, as apparent from the two m χ scenarios depicted in Fig. 5.In order to better understand this behavior, we introduce the "normalized" matrix elements for the W + W − γ channel, which is the only annihilation channel with non-vanishing diagonal and off-diagonal matrix elements in the wino/higgsino scenarios, where x ≡ E γ /m χ .These matrix elements are dimensionless and approach the LO fragmentation function of Eq. ( 21) asymptotically for sufficiently large m χ and small enough x.We plot them in Fig. 6.We observe that in the regime (1 − x) ≪ 1 and for scenarios with heavy

IJ
/ dx in the pure wino model, as obtained using our 'full SE' calculation (solid lines), with the 'endp SE' calculations of Refs.[84,85] (dashed lines), as functions of (1 − x) for three different values of m χ .Each color corresponds to the contribution of a single matrix element with indices IJ as introduced in Sec. 3. neutralinos, such as the considered case of m χ = 2.4 TeV, the diagonal matrix element with I = J = ⟨1 1⟩ is much larger than the corresponding matrix elements for all other IJ combinations.In this case the three relevant Sommerfeld matrix elements S IJ listed in Tab. 1 are comparable in order of magnitude.We see that, when they are combined with the relevant annihilation matrix elements, the last term of Eq. ( 22) in this scenario gives the dominant contribution in the (1 − x) ≪ 1 regime.The endpoint spectrum is thus dominated by virtual chargino annihilations.This is a consequence of the increasing influence of Sudakov logarithms for (1 − x) ≪ 1, which are particularly large for the terms of Eq. ( 22) involving charginos.
Besides the Sommerfeld effect, which is the focus of this work, a different kind of resummation is necessary at the endpoint region of the spectrum.Indeed, for very heavy neutralinos and E γ → m χ our 'full SE' computation suffers from large logarithmic enhancements that need to be incorporated (resummed) into the prediction.Such a resummation has already been completed for the MSSM in Ref. [62] at the next-to-leading logarithmic (NLL) accuracy, while in Refs.[61,84] and [80] an NLL-prime (NLL') accuracy8 was achieved for wino-and higgsino-like DM, respectively (see Ref. [87] for a short review).The NLL' calculation exhibits an accuracy of O(1%), while the accuracy of the NLL prediction is of order O(10%).Similar calculations for wino-like DM can be found in Ref. [88] at leading-logarithmic (LL) accuracy and Ref. [89] (NLL), and for higgsino-like DM in Ref. [90] (LL'), where all anomalous dimensions are given at 1-loop accuracy.
Indeed, we also verified that our results are consistent with the endpoint resummed calculation of Refs.[84,80,85], to which we refer as 'endp SE' in the following, within their kinematic regimes of validity, namely, m χ − E γ ∼ O(m W ). For these comparisons an additional step is required.While in the 'endp SE' calculations, each annihilation-matrix element [d( σv)/ dE γ ] IJ incorporates terms resummed to all orders in the electroweak coupling, annihilation matrix elements in the 'full SE' computation are computed at LO. Thus, for a meaningful comparison between the two calculations, a fixed-order expansion of each annihilation matrix element of the 'endp SE' calculation must be performed.The LO term of the fixed-order expansion of the IJ-th annihilation matrix element in the 'endp SE' calculation should then agree, up to power corrections of O(m 2 W /m 2 χ ) that are not included in our calculation, with the respective matrix element of the 'full SE' computation.Indeed, we find that the annihilation matrix elements associated with the χ 9 are exactly reproduced by the expanded results of Refs.[84,85] at the leading order in the power expansion assumed there, and provided that 1 − x ≪ 1.Note that, at LO, the γ + f f and γ + Zh final states can only occur when the initial state is composed of charginos, i.e. these channels only contribute to the diagonal I = J = ⟨1 1⟩ annihilation matrix element.In contrast, γ + W + W − final states can also result from neutralino annihilation, which receives contributions also from the IJ = ( 11) ( 11) and IJ = ( 11)⟨1 1⟩ combinations.
The agreement between our calculated annihilation matrix elements and those extracted from Refs.[84,85] is very good for masses above the TeV scale, but poor for lighter neutralinos.This is to be expected, as the 'endp SE' results are valid only up to power corrections of O(m 2 W /m 2 χ ) which are not negligible for smaller neutralino masses such as m χ = 150 GeV.Our calculations, instead, provide the most accurate and reliable picture for the gamma-ray spectrum from neutralino annihilations in that mass range.
In addition to the comparison with the 'endp SE' calculation, we compared our results for the I = J = ( 11) annihilation matrix element analytically with the fixed-order neutralino annihilation cross section into a W + W − γ final state in the pure wino and higgsino limits that has been provided in Ref. [81].We found exact agreement.

Conclusions
In this work, we have presented the first calculation of the continuum gamma-ray spectra resulting from neutralino annihilation including the Sommerfeld effect in the MSSM.The main novelty of our work is the systematic inclusion of all combinations of chargino-antichargino and neutralino pair annihilation processes into three-body final states that play a role in the calculation.The impact of the Sommerfeld effect is generically very strong due to the highly non-relativistic nature of the DM particles in the considered scenarios.For the sake of concreteness we focused our numerical discussion on the pure wino and pure higgsino limits of the MSSM.In the neutralino mass range of about 100 GeV to 1 TeV, we find qualitative differences compared to calculations of continuum spectra including only final-state radiation that are traditionally employed in gamma-ray searches for WIMP DM.For neutralinos heavier than about a TeV, the endpoint of the continuum requires the resummation of large Sudakov logarithms from soft/collinear electroweak radiation.Our work fills the gap between previous calculations focusing on the endpoint regime and separate ones for the low-energy photon regime where the widely used final-state radiation approximation is appropriate.
To ensure the correctness of our calculations we performed several stringent consistency checks.We verified that our results agree with older ones in the appropriate limits, we checked gauge invariance, unitary safety, and we re-derived the non-relativistic potential for the s-wave annihilation of neutralinos.We also showed that our results are consistent with existing results employing the collinear approximation.
On the technical side, our results will allow for a reassessment of the impact of internal bremsstrahlung by combining it with the Sommerfeld enhancement effect.This will pave the way for robust global fits, especially in reduced-parameter MSSM scenarios.Most importantly, though, our calculation will open the door to detailed phenomenological studies of the indirect detection of neutralinos using gamma-ray observations from both satellite and Cherenkov telescopes.In the light of improved energy resolutions and sensitivities of current and nextgeneration gamma-ray telescopes this is a very timely achievement.
In Eqs. ( 11)-( 16) the coupling matrices have been expressed in terms of the coefficients that are displayed below.In particular, the coefficients associated with massive vector boson interactions are given by v where c W = cos θ W and s W = sin θ W , and the N ij , Ũij and Ṽij denote elements of the unitary mixing matrices defined in Eq. ( 26).The coefficients that are associated with charged scalar bosons read s where c β ≡ cos β and s β ≡ sin β.Finally, the coefficients for the neutral scalar bosons are given by s s where s α ≡ sin α, c α ≡ cos α and α is the mixing angle of the CP-even Higgs doublet.

Figure 2 :
Figure 2: Prototypical ladder-like Feynman diagram contributing to the Sommerfeld effect.

Figure 6 :
Figure 6: Comparison of the normalized matrix elements dN W + W −