{'fraction_non_alphanumeric': 0.03684302427099703, 'fraction_numerical': 0.017932599965077702, 'mean_word_length': 4.379579184670298, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fermi LAT Search for Internal Bremsstrahlung Signatures from Dark Matter Annihilation

A commonly encountered obstacle in indirect searches for galactic dark matter is how to disentangle possible signals from astrophysical backgrounds. Given that such signals are most likely subdominant, the search for pronounced spectral features plays a key role for indirect detection experiments; monochromatic gamma-ray lines or similar features related to internal bremsstrahlung, in particular, provide smoking gun signatures. We perform a dedicated search for the latter in the data taken by the Fermi gamma-ray space telescope during its first 43 months. To this end, we use a new adaptive procedure to select optimal target regions that takes into account both standard and contracted dark matter profiles. The behaviour of our statistical method is tested by a subsampling analysis of the full sky data and found to reproduce the theoretical expectations very well. The limits on the dark matter annihilation cross-section that we derive are stronger than what can be obtained from the observation of dwarf galaxies and, at least for the model considered here, collider searches. While these limits are still not quite strong enough to probe annihilation rates expected for thermally produced dark matter, future prospects to do so are very good. In fact, we already find a weak indication, with a significance of 3.1\sigma (4.3\sigma) when (not) taking into account the look-elsewhere effect, for an internal bremsstrahlung-like signal that would correspond to a dark matter mass of \sim150 GeV; the same signal is also well fitted by a gamma-ray line at around 130 GeV. Although this would be a fascinating possibility, we caution that a much more dedicated analysis and additional data will be necessary to rule out or confirm this option.


Introduction
While there is little doubt about the very existence of a cold dark matter (DM) component in the present composition of the universe, contributing a substantial fraction of Ω χ = 0.229 ± 0.015 to its total energy budget [1], the DM nature remains unknown even 80 years after Zwicky first postulated a 'missing mass' in the Coma cluster [2]. At the moment, the leading hypothesis for a solution to this puzzle are thermally produced, weakly interacting massive particles (WIMPs) as expected to appear in many extensions to the standard model of particle physics (see, e.g., Refs. [3][4][5] for reviews or Ref. [6] for a recent pedagogical introduction)though one should keep in mind that the absence of any clear signs for new physics at the CERN LHC may render this theoretical prejudice less likely already in a few years from now [7]. If DM consists of WIMPs, it can be searched for both at colliders (with missing transverse energy as the main signature), in direct detection experiments looking for the scattering of WIMPs with the nuclei of terrestrial detectors, or indirectly through the observation of WIMP annihilation products in cosmic rays; the obvious advantage of indirect searches (see [8] for a recent status review) being that they allow to probe the nature of DM not only locally but on cosmological (or at least galactic) scales.
Gamma rays are for several reasons a particularly suitable detection channel for DM annihilation, not the least because they propagate essentially unhindered through the galaxy and thus directly point back to the sources. Here, we will focus on another important aspect, namely the possibility of sharp spectral features in the annihilation spectra that would allow a rather straightforward discrimination from astrophysical backgrounds. Indeed, it has early been pointed out [9] that the direct annihilation into photons would lead to the smoking gun signature of a gamma-ray line [10][11][12][13][14][15][16][17][18]. Equally pronounced spectral features near photon energies close to the kinematical endpoint at E γ = m χ arise due to internal bremsstrahlung, i.e. in the presence of three-body final states containing a photon [19][20][21][22][23][24], albeit at much higher rates because they are not loop-suppressed. Including such features in the analysis, rather than following the common approach of looking for featureless spectral templates, can significantly increase the sensitivity to DM signals [25,26] and might even be used to rather efficiently discriminate between different DM models (see e.g. Ref. [27]).
The Large Area Telescope (LAT) [28], the main instrument on board the Fermi Gammaray space telescope, has an unprecedented sensitivity to gamma rays from 30 MeV to above 300 GeV. This, together with its large field-of-view, makes it ideally suited for DM searches [29][30][31][32][33][34][35][36][37][38][39]; the non-observation of gamma-ray signals from dwarf galaxies, e.g., places the currently most stringent bounds on the annihilation rate of WIMPs with masses below around 700 GeV [30,31,39] (for higher masses, H.E.S.S. observations of the galactic center lead to stronger limits [40]). So far, only line-signals have been searched for in the LAT [36,37] or EGRET [41] data; the main purpose of this article is to extend these searches to other spectral endpoint features which, as outlined above, are arguably more relevant in most WIMP models. In order to look for such templates, we adopt a strategy that is very close in spirit to traditional line searches, as recently shown to be very promising [26], and use a new adaptive method to identify optimized target regions around the galactic center. As we will show, this leads to very competitive constraints on possible annihilation signals.
For Majorana (but also scalar) WIMPs, the pair of annihilating DM particles approximately forms a J = 0 state for the small relative velocities v ∼ 10 −3 expected in the galactic halo. In this commonly encountered situation, the annihilation rate into fermionic two-body final statesf f is (usually quite strongly) suppressed by m 2 f /m 2 χ and the next order result-with an additional photon in the final state-can be greatly enhanced by a factor of ∼ (α em /π)m 2 χ /m 2 f ; since this only happens for t-and u-channel annihilation mediated by the exchange of charged virtual particles, this process is also referred to as virtual internal bremsstrahlung (VIB, see Ref. [22] for an extensive discussion). Apart from the expected large rates, VIB also results in very pronounced spectral features close to the kinematic endpoint which resemble a slightly distorted line. Taken together, these two aspects make VIB in some sense the most promising DM signature to look for and motivate our choice of mainly focussing on related spectral distortions in the astrophysical background.
In order not to obscure our analysis by the potentially many parameters entering into the DM model, we will focus on a simple toy model that corresponds to the minimal extension to the standard model where we can expect strong VIB signals. While this model has been considered before and in its own right-in particular in connection with electroweak bremsstrahlung corrections to the annihilation rate [42][43][44][45][46][47][48], but also in other contexts [49,50]-let us stress that our analysis can straight-forwardly be applied to any other model with similar gamma-ray spectra from DM annihilation; in particular, the main features of this model are the same, for our purpose, as for neutralino DM in some phenomenologically very relevant regions of the supersymmetric parameter space.
The structure of this article is as follows. In Section 2, we introduce our toy model and comment on how it relates to the more commonly considered case of supersymmetric neutralino DM. We describe our method to search for pronounced spectral templates in the LAT data in Section 3 and also present our results there. In Section 4, we compare our new limits to existing limits from dwarf galaxies, expectations for thermally produced DM, collider constraints and limits from cosmic-ray anti-protons. We present our conclusions in Section 5. Finally, we provide some additional technical information about our method of selecting a target region optimized for the search of DM-related spectral features (Appendix A) and how a subsampling analysis of the full sky data can be used as a further test to confirm the reliability of our statistical method (Appendix B).

Toy model with large virtual internal bremsstrahlung
We will assume that the DM of the Universe is constituted by Majorana fermions χ, singlets under the Standard Model gauge group, which couple to the Standard Model via a Yukawa interaction with a scalar η that is not much heavier than the DM particle. The Lagrangian of the model reads: Here, L SM is the Standard Model Lagrangian. L χ and L η are the parts of the Lagrangian involving only the Majorana fermion χ and the scalar particle η, respectively, and are given by where D µ denotes the covariant derivative. Lastly, L int denotes the interaction terms of the new particles with Standard Model fields. We will consider in this paper three toy models where the DM particle only couples to the right-handed muons, tau leptons or bottom quarks, respectively, via a Yukawa interaction with the scalar η. We assume the latter to be an SU (2) L singlet in order to avoid constraints from electroweak precision measurements. The gauge quantum numbers of the intermediate scalar η are (1, 1) 1 for couplings with the muon or the tau (i.e. η is a SU (3) c and SU (2) L singlet with hypercharge Y = 1), and (3, 1) 1/3 for couplings with the bottom quark. Furthermore, to guarantee a coupling to just one generation of fermions we assign η a muon number L µ = −1, a tau number L τ = −1 or a beauty number B = −1, respectively. Then, the relevant part of the interaction Lagrangian reads with Ψ = µ, τ, b. Note that in principle additional couplings of the form H † Hη † η and (η † η) 2 are allowed (where H denotes the Higgs doublet). We will neglect them throughout this work since they do not directly influence the gamma-ray signature we are interested in. In these scenarios, DM particles can annihilate into two fermions with a velocityweighted annihilation cross-section which can be decomposed into an s-wave and a p-wave contribution. The s-wave contribution reads in lowest order of the relative center-of-mass velocity v [51,52] (σv) s-wave 2-body = where µ ≡ (m η /m χ ) 2 parametrizes the mass splitting between the DM particle χ and the t-channel mediator η, and the color factor N c is one for muons and taus and three for bottom quarks. The p-wave contribution at lowest order in v is [50] (σv) p-wave It is important to note that the s-wave contribution to the velocity-weighted annihilation cross-section of Majorana fermions is helicity suppressed, by the mass squared of the daughter fermion, whereas the p-wave contribution is suppressed by the velocity squared of the galactic DM particles today, typically v ∼ 10 −3 . Therefore, the 2 → 2 annihilation cross-section is fairly small, and higher order corrections could play an important role.
Indeed, it was shown in Refs. [53,54] that the associated emission of a vector boson lifts the helicity suppression in the s-wave contribution to the annihilation cross-section. This process was later dubbed virtual internal bremsstrahlung (VIB), which together with photons from final-state radiation (FSR) constitute the full internal bremsstrahlung (IB) spectrum [22]; the corresponding Feynman diagrams are shown in Fig. 1. The explicit expression for the annihilation cross-section into two massless fermions and one VIB photon is [45,55]: which is usually non-negligible and can in certain instances, in particular for small values of µ, be considerably larger than the 2 → 2 annihilation cross-sections Eqns. (2.4,2.5); Q f denotes the electric charge of f and η in units of |e|. We emphasize that throughout this paper "3-body process" refers to the VIB process only, whereas "2-body process" refers to the helicity-suppressed tree-level process χχ → ff plus the FSR photons. When explicitly referring to (σv) 2-body in the following, we will therefore multiply Eqs. (2.4,2.5) by a factor of 1 + dx dN FSR /dx , where [20] Furthermore, the energy spectrum of gamma rays produced in the 2 → 3 process has a very peculiar shape that allows for an efficient search for gamma rays from internal bremsstrahlung. Namely, the differential three-body cross-section, as function of the VIB photon energy x ≡ E/m χ , is given by [22] v To illustrate the peculiar features of VIB, the energy spectrum of gamma rays that is produced per annihilation in our toy model is shown in Fig. 2 for the three different final state fermion flavours; for definiteness we assume m χ = 200 GeV and a relatively small  Figure 1. Feynman diagrams of the processes that contribute in leading order to the three-body annihilation cross-section and produce internal bremsstrahlung. The first diagram very roughly corresponds to VIB, the second and third to FSR (but note that these contributions can be properly defined and separated in a gauge-invariant way [22]). mass-splitting of µ = 1.1. The spectra of secondary photons that stem from the subsequent decay or fragmentation of the produced fermions are derived using Pythia 6.4.19 [56]. Note that in case of bottom-quark final states we also take into account the production of VIB gluons following Refs. [48,57]. 1 For two-body annihilation, we cross-checked our results with the analytical fits from Ref. [58,59] and find very good agreement. From Fig. 2 it is clear that for small enough mass-splittings the gamma-ray spectrum at high energies is completely dominated by VIB photons, which show up as a pronounced peak at energies close to the dark matter mass. Secondary photons and FSR only become relevant at lower energies, or for larger values of µ. In our spectral analysis of galactic center fluxes presented in Section 3, we will entirely concentrate on the spectral VIB feature and neglect the featureless secondary photons. We will consider the range 1 < µ 2, because the VIB feature is most important in the nearly degenerate case. In this range, the shape of the VIB spectrum is almost independent of µ (it becomes slightly wider for larger µ), but its normalization can vary rather strongly: for µ = 1.1 (µ = 2.0), the rate is already suppressed by a factor of 0.55 (0.05) with respect to the exactly degenerate µ = 1 case; for large µ, the rate scales as ∝ µ −4 (whereas the two-body annihilation rate scales like ∝ µ −2 ). For comparison with our main results, we will also derive limits from dwarf galaxy observations (see Section 4.1); in this case we will take into account both VIB and secondary photons.

Connection to the MSSM
Before continuing, let us briefly mention the connection between our toy model and the much more often studied case of supersymmetry. The minimal supersymmetric extension to the standard model (MSSM) is extremely well motivated from a particle physics point of viewleading, in particular, to a unification of gauge couplings and strongly mitigating fine-tuning issues in the Higgs sector-and the stability of the lightest supersymmetric particle (LSP) is guaranteed by the conservation of R-parity; if it is neutral and weakly interacting, the LSP thus makes for an ideal DM candidate (for a comprehensive and pedagogical primer to supersymmetry and the MSSM see e.g. Ref. [61]).
In most cases, the lightest neutralino is the LSP, and thus a prime candidate for WIMP DM [3]. It is a linear combination of the superpartners of the neutral components of the  U (1) × SU (2) gauge as well as Higgs fields, and thus a Majorana fermion just like the DM particle in our toy model. As pointed out above, the annihilation into fermion-antifermion pairsf f is therefore helicity suppressed in the limit of small velocities; this helicity suppression can be lifted if an additional photon is present in the final state and annihilation happens via the t-channel exchange of a charged particle. In the case of supersymmetry, this can only be achieved through the corresponding left-and right-handed sfermionsf L andf R which, in the limit of vanishing m f , couple to the neutralino and fermions as where as usual f R/L ≡ 1 2 (1 ± γ 5 )f . Compared to Eq. 2.3, the sfermions thus play exactly the same role as η and the main difference to our toy model is that i) there are two relevant scalars for each fermion final state and that ii) the interaction strength y (R,L) is no longer a free parameter but uniquely defined by gauge symmetry, and of course the composition of the neutralino (see e.g. Ref. [62]): where g is the usual SU (2) coupling and θ W the Weinberg angle. In Eq. (2.11), the upper (lower) signs apply if the third component of the weak isospin for f is given by While the couplings are fixed, the relative importance of thef f γ final state to the total annihilation rate in the v → 0 limit strongly depends, as we have seen, on the mass difference µf ≡ (mf /m χ ) 2 between neutralino and sfermion. In a phenomenological description of the MSSM, i.e. without referring to a specific mechanism for supersymmetry breaking, all µf can be treated as essentially free parameters; choosing one sfermion to be much more degenerate in mass with the neutralino than all the others thus effectively, for the purpose of our discussion, reduces the general MSSM case to our toy model. 2 In fact, such a situation is phenomenologically very relevant already in the most minimal supersymmetric setup, the constrained MSSM (CMSSM), and appears both in theτ -and (for neutralino masses larger than what we are interested in here) in thet-coannihilation region-which apart from the focus point, funnel and bulk region are the only regions in the CMSSM parameter space where the neutralino acquires the correct thermal relic density (see e.g. Ref. [63] for a discussion; note, however, that the bulk region now is already excluded by LHC data). Correspondingly, it was found in Ref. [22] that VIB can indeed greatly dominate the total photon spectrum at high energies in the MSSM in general, but in particular in the coannihilation regions of the CMSSM. If only the stau is degenerate in mass with the neutralino (as is the case for example in benchmark model BM2 of that reference), the shape of the resulting gamma-ray spectrum at high energies is almost indistinguishable from the τ + τ − γ case displayed in Fig. 2; the same holds for the more common situation that all leptons are rather degenerate with the neutralino (as in benchmark model BM3; small shape differences in the high energy spectra can be attributed to slightly different mass splittings). BM3 lies with a neutralino mass of m χ = 233.3 GeV in the reach of Fermi LAT, and we will comment on this particular benchmark point in light of our results below.

Fermi LAT search for VIB signatures at the galactic center
The VIB gamma-ray signal produced in our toy model is sharply peaked at energies close to the DM mass, very much like gamma-ray lines that are produced in the two-body annihilation into γγ. For this reason, our analysis will closely follow the spirit of traditional gamma-ray line searches [36,37,41], as some of us recently proposed in Ref. [26]. The practical advantage of DM signals that are very concentrated in the energy spectrum is that searches for these signals can be restricted to relatively small energy windows, being only a factor of a few larger than the energy resolution of the instrument. For such small energy ranges, it is reasonable (and in fact a posteriori justified by our results) to approximate astrophyscial background fluxes by a simple power-law; its normalization and spectral index are obtained directly from a fit to the data. Since a detailed understanding of the background sources is not necessary in this case, it is possible to choose even complex (but very promising) regions like the Galactic center as target in the search for a DM signal, as we will do in our present analysis.
We will here concentrate on DM masses in the range 40 GeV < m χ < 300 GeV. The lower end of the mass range is motivated by the LEP constraint on the mass of charged exotic particles, which reads m η > 40 GeV (see Section 4.3); since we are mainly interested in the degenerate scenario where m χ ≈ m η , this already excludes much lighter dark matter particles. For DM masses above 300 GeV, on the other hand, the spectral feature would be outside of the nominal energy range of the Fermi LAT. As discussed above, we concentrate here on mass splittings in the range 1 < µ 2, for which the spectral shape of the VIB signal is practically independent of µ.

Dark matter signal from the galactic halo
The gamma-ray flux from DM annihilation in the galactic DM halo is given by a line-of-sight integral over the DM density squared, Here, m χ is the DM mass, σv the total DM annihilation cross-section averaged along the line of sight, dN/dE the energy spectrum of produced gamma rays, and ξ denotes the angle to the Galactic center. The coordinate s ≥ 0 runs along the line of sight, and the distance to the Galactic center r is given by r(s, ξ) = (r 0 − s cos ξ) 2 + (s sin ξ) 2 , where r 0 = 8.5 kpc denotes the distance between Sun and the Galactic center. We will consider the following generalized Navarro-Frenk-White (NFW [36,64]) profile normalized to the fiducial value ρ χ = 0.4 GeV cm −3 at Sun's position [65,66] and with a scaling radius of r s = 20 kpc. In case of an inner slope of α = 1 this reproduces the standard NFW profile. The possible impact of adiabatic contraction [67][68][69][70] can be studied in an effective way by allowing for larger inner slopes of the profile. We will concentrate on the range 1 < α 1.4, which is still compatible with microlensing and dynamical observations [71] (traditional adiabatic contraction following Ref. [67] would give rise to even larger values of α, see e.g. Ref. [71]). The modified isothermal and Einasto profiles [72][73][74] are expected to give comparable results to the standard NFW profile in searches for line-like features [37].

Event and target region selection
The gamma-ray data measured by the Fermi LAT is publicly available [75]. The events that enter our main analysis are taken from the P7CLEAN V6 event class. We consider 43 months of data (from 4 Aug 2008 to 6 Feb 2012), and select front-and back-converted events with energies in the range 1-300 GeV. We apply a zenith angle cut of θ < 100 • in order to avoid contamination with photons from the earth albedo, as well as the quality cut filter DATA QUAL==1 (all event selection is done using the 06/10/2011 version of ScienceTools v9r23p1. The target region for our spectral analysis is chosen by using a data-driven adaptive procedure with the aim of maximizing the expected signal-to-noise ratio. We stress that such an approach is extremely important when looking for spectral features at the statistical limit of the detector, and that an inefficiently chosen target region can easily wash out or hide a potential signal. Our choice of the optimal target region depends on the adopted DM profile; in order to select it, we estimate the expected spatial distribution of background noise in our search for spectral features above 40 GeV by considering the actually measured events below 40 GeV. The spatial distribution of signal photons, on the other hand, just follows from Eq. (3.1). All details of the method are given in Appendix A. We adopt four reference values for the inner slope of the DM profile, α = 1.0, 1.1, 1.2 and 1.4, for which we obtain the target regions that are shown in Fig. 3 as solid black lines. In this plot, the colors encode the expected signal-to-background ratio in different regions of the sky, normalized to one for the pixel where this ratio is maximal (note that the actual value of this quantity is a factor of 1.9 (3.9, 31) larger for Reg2 (Reg3, Reg4) than for Reg1). In case of a standard NFW profile with α = 1.0, the target region includes besides the galactic center also regions at higher and lower latitudes up to |b| 70 • ; for steeper profiles the optimal target regions shrink drastically to regions closer to the galactic center. The galactic disc is strongly disfavoured in all cases. Southern regions are somewhat preferred, since the diffuse gamma-ray emission from our galaxy is not perfectly north/south symmetric. From these four regions we extract the measured spatially integrated gamma-ray energy spectrum for our subsequent analysis.

Spectral analysis
The search for VIB signatures is done by using the sliding energy window technique discussed e.g. in Refs. [26,36,37,41]: we consider for each DM mass m χ in the range 40 GeV < m χ < 300 GeV a small energy window that is approximately centered on m χ , and hence on the position of the expected VIB feature. More precisely, we use the energy range E = m χ ǫ −0.7 . . . min[m χ ǫ 0.3 , 300 GeV], where the size of the window ǫ varies between ǫ ≃ 1.8 for m χ = 40 GeV and ǫ ≃ 7 for m χ = 300 GeV. The size of the window is identical to the values used in Ref. [37], where it was found to lead to reasonable background fits in context of gamma-ray line searches. The position of the window is not exactly symmetric around m χ , but slightly shifted towards lower energies as it was suggested for VIB features in Ref. [26] in order to increase the sensitivity. We emphasize again that secondary photons, as they come from the decay or fragmentation of the fermions, are neglected in our spectral analysis, because they become relevant only outside of the energy window that we consider here (at least for our toy model).
For each given mass m χ , and within the adopted small energy window, we now fit the gamma-ray spectra measured in the different target regions of Fig. 3 with a simple threeparameter model: The astrophysical background fluxes are approximated by a power law with a free spectral index (1) and normalization (2); the DM VIB signal has only a free normalization (3), whereas its mass and the mass-splitting (which we set to µ = 1.1 in most of the analysis) remains fixed during the fit. For physical reasons we require the normalization of the VIB signal to be positive.
Technically, we perform a binned analysis of the gamma-ray spectrum measured in the different target regions. To this end, we distribute the corresponding measured gamma-ray events in a very large number of energy bins (200 per decade). Since the size of the adopted energy bins is much smaller than the energy resolution of the Fermi LAT, our analysis is essentially identical to an unbinned analysis of the energy spectrum. Each energy bin j then contains a number c j of events. The number of expected events µ j are obtained by convolving our above three-parameter model with the energy dispersion and the exposure of the LAT; the resulting µ j can then be fitted to the observed counts c j by maximizing the likelihood function L = Π j P (c j |µ j ) with respect to the three model parameters. Here, P (c|µ) is the Poisson probability to observe c events when µ are expected. Note that the functional form of the energy dispersion is directly inferred from the IRF of the P7CLEAN V6 event class and correctly averaged over impact angles and front-and back-converted events, using our own software. Exposure maps are derived using the ScienceTools v9r23p1.
Limits on or the significance of a DM VIB contribution can then be derived by using the profile likelihood method [76]. A one-sided 95% C.L. upper limit on the annihilation cross-section is obtained by increasing the DM signal normalization from its best-fit value until −2 ln L increased by a value of 2.71 (while refitting the background parameters). The significance of a signal, on the other hand, is derived from the test statistics (TS) where L best-fit is the likelihood of the best-fit model, and L null the likelihood of the null hypothesis (the absence of a DM signal; the null model has only two free parameters). In absence of a signal, one expects that the TS follows some χ 2 distribution. More precisely, since the normalization of the DM signal is bounded to be positive, the TS should follow a 0.5χ 2 k=0 + 0.5χ 2 k=1 distribution (see e.g. Ref. [  degree of freedom, respectively. 3 This theoretical expectation is indeed very well confirmed by a subsampling analysis of the data as we discuss in Appendix B. Hence, if the test statistics is measured to be T S for a certain DM mass in a single trial, this would correspond to slightly more than √ T Sσ significance. However, since in the present analysis we effectively perform many statistically independent trials when scanning through m χ and analyzing different target regions, the probability to find just by chance a statistical fluctuation that mimics a signal is increased; this is known as the look-elsewhere effect (LEE). In our case, we approximate the distribution of maximal TS values from which the significance of our signature is calculated by 4 × 4 = 16 trials over a χ 2 k=2 distribution. Four trials over a χ 2 k=2 distribution come from the scan over m χ (see Appendix B or, for a general discussion, Ref. [77]); the remaining four trials are associated with the four target regions. In practice, the significance of the observed signature is then found by solving P (χ 2 k=2 < T S) #trials = P (χ 2 k=1 < σ 2 ) for σ, where P (χ 2 k < x) denotes the probability to observe a value smaller than x when drawing from a χ 2 k distribution.

Results
In our analysis, no VIB signal with a significance of at least 5σ were found. Instead, we show in Fig. 4 the upper limits at 95% CL on the VIB three-body cross-section, for the different values of the inner slope α that correspond to the target regions in Fig. 3. We assumed a mass-splitting of µ = 1.1 for definiteness; limits for other values of µ are very similar and will be discussed below (see e.g. Fig. 6). As can be seen from the plot, our limits are always stronger than the 'thermal cross-section' of 3 × 10 −26 cm 3 s −1 that is often quoted for comparison; in the case of contracted profiles with α = 1.4 they can even reach down to values of 10 −28 cm 3 s −1 for DM masses m χ 100 GeV. As we will discuss below in Section 4.1, our limits are much stronger than what can be obtained from e.g. dwarf galaxy observations. For comparison, the gray cross in Fig. 4 shows the CMSSM benchmark point BM3 [22], which lies in the coannihilation region and was already discussed above. This benchmark point still remains unconstrained by our limits; its rather small cross-section is closely related to the requirement that the neutralino is a thermal relic, as we will discuss in Section 4.2 below.
In the left panel of Fig. 5, we show the significance for a VIB-like spectrum as function of m χ , assuming that µ = 1.1. The different lines correspond to the different target regions. The significance is shown in terms of the TS value that was discussed above. We find a possible signal candidate at a DM mass of m χ ≈ 150 GeV. The indication for a signal is largest for the target region Reg2, which corresponds to α = 1.1, and has a nominal significance of √ T S = 4.3σ. Taking into account the LEE as discussed above, the significance is 3.1σ. The corresponding fit to the data is shown in the right panel of Fig. 5; the spectral feature in the measured flux can be easily recognized by eye. A similar preference for a signal, although with less significance, appears also in the other regions Reg1, Reg3 and Reg4 (note that the fluctuations around 50 GeV are completely within the statistical expectations). TS values of zero indicate that for these values of m χ the data would be best fitted with an unphysical negative signal normalization; in this case, the likelihood of the model with DM contribution becomes identical to that of the null model because we enforced a non-negative signal normalization in our fits.
We have performed several tests to exclude the tempting DM interpretation of this signature, none of which has succeeded so far: By masking out different halves of the signal region of Reg2, for example, we find that the signal independently appears in the north, south, east and west parts of Reg2 (though with a large scatter in the significances), as expected from a DM signal. When shifting the target region away from its position by about 10-20 • , on the other hand, the signal disappears completely. This makes a purely instrumental effect, which would likely also appear in other sky regions, less likely. To test our statistical method, we performed a subsampling analysis of the Fermi LAT data in the galactic anticenter region in order to obtain an empirical understanding of the statistical behaviour of our TS (see Appendix B for details). We find that-in absence of a signal-our TS follows very well a χ 2 distribution, as theoretically expected. The signature appears also in the P7ULTRACLEAN V6 and P7SOURCE V6 event classes, and when considering back-and front-converted events separately; furthermore, as expected for a real signal, its significance has grown with time (i.e. we find smaller T S values when considering instead older and older data sets, though with some scatter around the linear trend). The uncertainty in the effective area as relevant for searches for line-like features is about 2% [78], which is much smaller than the fractional contribution of the observed signature to the flux.
Despite these encouraging facts, we call for caution when interpreting this signature as due to DM annihilation. First of all, with a significance of 3.1σ (when taking the LEE into account), it could still simply be an upward fluctuation at the right place. Alternatively, the observed signal could be real but due to a yet unidentified astrophysical process, like e.g. the inverse Compton scattering of an extremely hard electron component on stellar light (see e.g. Ref. [79] for a related discussion of ICS in the context of the Fermi bubbles)though in general, as already stressed several times, it is quite difficult to explain this kind of spectral features with astrophysical processes. In this context, it is also worth to mention again that our analysis crucially depends on the commonly adopted assumption that the background locally follows a power-law, i.e. within each energy window that we consider. In principle, it might thus make for an interesting follow-up study to perform a signal fit on a more complicated background which contains, e.g., a break in the spectral index that could be confused with a signal. We do not expect, however, that such an analysis leads to qualitatively different results than presented here because the data itself tells us that the single power-law assumption works very well (see Appendix B), and because of the sharpness of the observed signature. Even more importantly, we always find a spectral index for the background contribution that is consistent with −2.6; this value is expected for the production of gamma rays from the collision of cosmic rays with the interstellar medium and thus extremely well motivated from astrophysics (see e.g. Ref. [80]). Lastly, our analysis relies entirely on the publicly available data, which makes it impossible to take into account all known instrumental effects. However, we strongly believe that our actual statistical analysis is sound, and not significantly affected by the obvious systematics of the LAT.
Finally, we note that the best-fit three-body cross-section for a VIB signal in Reg2 is σv = (6.2 ± 1.5 +0.9 −1.4 ) × 10 −27 cm 3 s −1 (assuming α = 1.1), the best-fit mass is in the range m χ = 149 ± 4 +8 −15 GeV; the errors are respectively statistical and systematical. 4 This is quite a bit larger than what is actually expected for VIB from thermally produced DM (see Section 4.2 for a discussion), making a straightforward interpretation of the signal in terms of a VIB signal somewhat less appealing. However, let us point out that the signature can also well be fitted by a pure gamma-ray line as produced in χχ → γγ, with a dark matter mass of around m χ ≃ 130 GeV and an annihilation cross-section σv ∼ 10 −27 cm 3 s −1 (see also Ref. [81] for a dedicated analysis). Finally, we note that the more commonly adopted annihilation spectra from e.g. χχ → µ + µ − or χχ →bb annihilation are much too flat to fit the data and thus cannot be used to explain the signal.   Figure 7. Same as Fig. 6, but for annihilation into τ + τ − γ.

Comparison with other constraints
The limits on the cross-section σv 3-body as shown in Fig. 4 were derived from VIB photons only, according to the spectrum given in Eq. (2.8). We will now discuss these constraints in light of other complementary probes, namely limits from gamma-ray observations of dwarf galaxies, limits that can be derived from the thermal production of DM and collider limits; we finally stress that small values of µ, for the annihilation intobbγ final states, can also very efficiently be probed by both cosmic ray antiprotons and direct DM detection experiments.  Figure 8. Same as Fig. 6, but for annihilation intobbγ. Note that we includebbg final states when calculating the dwarf limits [48,57].

Limits from dwarf galaxies
Important limits on the DM annihilation cross-section derive from gamma-ray observations of dwarf galaxies [30,31,39]. To compare our above results in Fig. 4 with the constraints that typically can be obtained from these targets, we will make use of the flux limit presented in Ref. [39]. 5 It reads, in terms of the total annihilation cross-section, σv < 8π m 2 χ /N tot γ · 5.0 × 10 −30 cm 3 s −1 GeV −2 . Here, N tot γ denotes the number of photons that are produced per annihilation in the energy range 1-100 GeV. When calculating N tot γ , we fully take into account all photons from three-body as well as the common two-body final states according to Eqs. (2.4) and (2.6); 6 the relative importance of two-and three-body final states depends on µ and the mass of the final state fermions.
We plot the resulting limits in Figs. 6, 7 and 8 as green lines for different values of the mass splitting µ and for the three different final state fermion flavours; in order to allow a simple comparison with our result in Fig. 4, we choose to present them in terms of the threebody annihilation cross-section (by rescaling them by a factor of σv 3-body /[ σv 2-body + σv 3-body ]). In most cases the limits depend relatively strongly on the mass splitting µ (Figs. 7 and 8); only for annihilation into µ + µ − the limits remain practically independent of µ, since the fluxes are always dominated by VIB photons (cf. Figs. 2 and 6). In contrast to that, our limits from the spectral search for VIB features in the galactic center, as shown by the black lines, exhibit only a very weak dependence on the mass splitting parameter µ for values of µ ≈ 1.1-2, since the spectral shape of the VIB radiation changes only mildly in this range.
As long as the mass splitting µ remains small, and for the leptonic final states, our spectral search in the galactic center fluxes for VIB features leads to much stronger constraints on our toy model than dwarf galaxy observations, see Figs. 6 and 7. For a large enough mass splitting, or in case of colored final state fermions like bottom-quarks, however, the 3-body annihilation into VIB photons becomes less relevant and the dwarf limits start to overtake our VIB limits (cf. Fig. 8). Note that for colored final states the VIB of gluons becomes important [48,57] and actually dominates the dwarf limits (see also Section 4.4).

Thermal production
In our toy model, self-annihilation in the early Universe is usually dominated by the p-wave process χχ →f f , since the averaged velocity is v ∼ O(1), see Eq. (2.5). Any embedding of this model in a full theory (like the MSSM, see Section 2.2) would likely increase the total annihilation rate. A loose but reasonable upper limit on the coupling constant y can hence be derived by requiring that the relic density which follows from the t-channel process χχ →f f alone does not already undershoot the measured DM density. From this, we could estimate an upper limit on the partial annihilation cross-section for χχ →f f γ as follows: At the freezeout temperature T f ∼ m χ /20, the p-wave annihilation channel χχ →f f is mildly suppressed by a velocity factor of v 2 = 6T f /m χ ∼ 0.3; the corresponding annihilation cross-section should be equal to or smaller than the thermal cross-section σv th ≃ 3 × 10 −26 cm 3 s −1 , in order to not undershoot the observed DM relic density. However, the s-wave VIB process χχ →f f γ is suppressed by a factor ∼ α em /π (for final state fermions with charge one), independently of the temperature T or the velocity v. Taking these pieces of information together, one finds an upper limit of roughly αem 0.3π σv th ∼ 2 × 10 −28 cm 3 s −1 for the threebody annihilation cross-section σv 3-body . Hence, if we want our toy model to be compatible with the observed relic density (assuming thermal DM production and a standard thermal history of the universe), this suggests an annihilation cross-section σv 3-body that is about two orders of magnitude smaller than the common thermal cross-section. Note, however, that e.g. destructive interference with s-channel diagrams, which do not appear in our toy model, could spoil this argument and thus in principle allow for larger rates.
Neglecting the coannihilation of χ and η, the relic density due to two-body annihilation χχ →f f would be approximately given by (see e.g. Ref. [46]) For comparison, we performed a full calculation with MicrOMEGAs [83]. We find that our order-of-magnitude estimate is only approximately valid if µ 1.3 in case of leptonic final states, and if µ 1.8 in case of bottom quark final states (see also Ref. [82]). For smaller mass splittings, coannihilations play a significant role. Note furthermore that even a vanishingly small coupling y may in fact not sufficiently suppress the total effective annihilation rate, for low dark matter masses m χ , to reproduce the observed relic density simply because the gauge-boson mediated annihilation of ηη pairs will at some point start to dominate. This effect is most pronounced for small mass differences µ and can clearly be seen in Figs. 6-8 (see also Ref. [84]). Furthermore, the coupling to the Higgs can become relevant at masses around 50 GeV. Using our full numerical results, and requiring now that Ω χ h 2 0.1, we obtain upper limits on the cross-section as shown in Figs. 6-8 by the blue lines for different values of µ. These limits largely agree with the rough estimates given in the preceding paragraph as long as the mass-splitting is µ ∼ 1.5. For larger mass splittings, or for smaller splittings where coannihilation would further reduce the relic density, the limits become stronger. From Figs. 6-8 it is apparent that in case of a regular NFW profile with (α = 1) our above galactic center limits are still around one (two) orders of magnitude away from values of the three-body cross-section that are consistent with a thermal relic in case of lepton (bottom quark) final states. However, in case of compressed profiles with α > 1, our limits become stronger, as shown in Fig. 4, and for inner slopes of the DM halo α ≈ 1.2-1.4, they can become sensitive enough to probe the required cross-sections of ∼ 10 −28 cm 3 s −1 and below.
Further enhancements of the signal, without violating our relic density constraints, can appear due to the gravitational clustering of DM as predicted by N -body simulations [74,85]. Moreover, the effect of coannihilation could actually be inversed and increase the relic density in some cases [86], which again would allow larger annihilation rates while still being compatible with the relic density. Finally, we note that the non-thermal production of dark matter or a non-standard freeze-out history would invalidate our relic density limits and admit larger annihilation cross-sections today (see e.g. Refs. [87][88][89]).

Limits on light charged or colored scalar particles
Scenarios with light scalars that couple to the Standard Model quarks and leptons are constrained by direct searches at colliders as well as by indirect searches through their effect on the quantum corrections to the electroweak precision measurements.
Collider searches. The most severe constraint on scenarios with light charged particles follows from the non-observation of an excess in the invisible decay width of the Z boson at LEP, ∆Γ inv < 2.0 MeV [90], which excludes the existence of exotic charged scalar particles with mass below 40 GeV [91]. Since we are mostly interested in the degenerate case where m χ ≈ m η , we only considered DM masses with m χ > 40 GeV in this work.
Additional constraints on the mass of charged scalar particles can be derived from the non-observation at LEP of an excess over the Standard Model expectations of dilepton events with missing energy, generated by the production of charged scalar particles that decay into a lepton and an invisible particle (in our framework, the DM particle). The search strategy relies on the identification of at least one lepton event as well as on satisfying the requirements on isolation and transverse momentum. Hence, the efficiency of the search decreases drastically when the DM particle and the intermediate scalar are more and more degenerate in mass, which is exactly the region we are interested in. The best constraints on such a scenario were published by the OPAL collaboration using 680 pb −1 of e + e − collisions at center-of-mass energies between 192 GeV and 209 GeV; they read [92] • m η ≥ 94.0 GeV at the 95% CL for couplings to muons, assuming m η − m χ > 4 GeV and 100% branching ratio for η → µχ • m η ≥ 89.8 GeV at the 95% CL for coupling to taus, assuming m η − m χ > 8 GeV and 100% branching ratio for η → τ χ.
For µ ≈ 1.1, these limits are always satisfied by our toy model. However, much larger mass splittings are already partially ruled if the dark matter mass is below ∼ 90 GeV. Searches for light colored scalar particles have been undertaken both at LEP, Tevatron and at the LHC. The non-observation of an excess over the Standard Model background of dijet events with missing energy translates into the following 95% CL limits on the mass of colored scalars, assuming 100% branching ratio for η → bχ (see also Fig. 4 in Ref. [57] for an updated summary): • m η ≥ 875 GeV, assuming m η − m χ > 130 GeV (ATLAS [93]) • m η ≥ 76 GeV, assuming m η − m χ > 7 GeV (DELPHI [94]) • m η ≥ 89 GeV, assuming m η − m χ > 8 GeV (ALEPH [95]).
Again, as long as µ ≈ 1.1, our toy model is not affected by these limits; only much larger mass splittings are already partially ruled out.
Oblique parameters. The new scalar particles carry charges under the electroweak group and can in principle modify the oblique parameters S, T and U [96,97]. However, in the scenario where the scalar particle is an SU(2) L singlet, none of the three oblique parameters receives any exotic contribution.
Muon g − 2. The interaction which leads to the VIB process χχ → µ + µ − γ leads, upon closing the DM fermion loop, to a contribution to the muon g − 2. In the toy model presented in this work, the contribution from DM particles to the muon g − 2 reads (see e.g. [98]) where µ ≡ (m η /m χ ) 2 denotes the usual mass splitting, which enters the function This contribution should be compared to the 3.2σ deviation between the Standard Model prediction and the experimental measurement [99], which is of opposite sign. Therefore, in our toy model the discrepancy between the theoretical prediction of the muon g−2 with the experimental measurement is larger than in the Standard Model. If we interpreted the g − 2 anomaly as a statistical fluctuation, the total theoretical prediction should still not deviate more than 5σ from the experimental value, which implies in the limit µ → 1 an upper bound on the coupling of y 1.7(m χ /100 GeV). The corresponding upper limit on the three-body annihilation cross-section reads (4.5) Comparing this to our VIB limits shown in Fig. 4, we find that-within our toy-model and the above assumptions-the g − 2 constraints are typically weaker than what we get from the spectral analysis of gamma-ray fluxes. Alternatively, additional particles could exist that generate a positive contribution to the muon g − 2 which compensates the negative contribution from the DM particles, thereby bringing the theoretical prediction closer to the experimental value.

Anti-proton observations and direct WIMP searches
In our toy model, the presence of a light colored scalar η opens up the gluon-VIB channel χχ →f f g; in general, this process has a ∼ 100 times larger cross-section than χχ →f f γ, because instead of α em the strong coupling α s enters the diagrams [48,57,100,101]. Hence, the limits on χχ →bbγ that can be obtained by constraining the corresponding process χχ →bbg with cosmic-ray anti-proton observations are quite significant for small values of µ; depending on the details of cosmic-ray propagation, they can be comparable to or even stronger than our gamma-ray limits from the VIB search [48,57]. Furthermore, χχ →bbg contributes to the gamma-ray dwarf limits in Fig. 8 [57]. The leptophilic models discussed in this paper also give rise to an antiproton flux through the annihilation channel χχ →f f Z which, when kinematically open, has a cross-section comparable to χχ →f f γ [42][43][44][45][46][47][48]. For these models, our gamma-ray limits are one order of magnitude stronger than the antiproton limits in the mass range m χ ∼ 100-300 GeV.
Lastly, for our bottom quark scenario, limits from direct DM searches may become very relevant, as the scattering between χ and the detector nuclei via an intermediate η becomes resonant for small values of µ [102].

Conclusions
One of the main challenges that are encountered when searching for signatures from DM annihilation in the cosmic gamma-ray flux is the discrimination between a possible signal and the astrophysical background. Very pronounced spectral features, like gamma-ray lines or photons from internal bremsstrahlung, play a key role in this context because it would be difficult to attribute them to astrophysical processes. Previous searches mostly concentrated on gamma-ray lines; here, we instead analysed the effects of internal bremsstrahlung which is intrinsically more promising due to the larger rates that are expected.
To this end, we defined a simple toy model (which can be considered as a subset of the MSSM particle content) with the important feature to generate an intense gamma-ray signature from virtual internal bremsstrahlung. We performed a dedicated search for this signature in 43 months of data from the Fermi Large Area Telescope, building on standard methods for gamma-ray line searches. In order to determine the optimal target regions for our spectral analysis, we introduced a new adaptive method that takes into account different conventional and cuspy DM halo profiles, see Fig. 3. We believe that this method will turn out to be useful even in other contexts, essentially whenever spectral features are being looked for at the statistical limits of the detector. For our toy model, in case of leptonic final states, we find upper limits on the annihilation cross-section that are stronger than what can be obtained from dwarf galaxies or collider searches, see Figs. 4, 6, 7 and 8. Our limits are still about an order of magnitude too weak to constrain annihilation cross-sections naively expected for a thermal relic (assuming a standard NFW profile for the DM density and no enhancement of the annihilation rate due to substructures). However, future prospects to do so are quite good [26].
We also find a weak indication for an internal bremsstrahlung signal, see Fig. 5, corresponding to a DM mass of m χ = 149 ± 4 +8 −15 GeV and an annihilation rate of σv χχ→f f γ = (6.2 ± 1.5 +0.9 −1.4 ) × 10 −27 cm 3 s −1 (we note here that the same signal can also be fitted with a conventional gamma-ray line at around 130 GeV and with a cross-section of about σv χχ→γγ ∼ 10 −27 cm 3 s −1 , and refer to Ref. [81] for a more detailed discussion). After taking into account the look-elsewhere effect, the signal significance is about 3.1σ (without, it is 4.3σ). We have performed several statistical tests of our method and deem that a purely instrumental effect would be a very unlikely explanation for this signal. As stressed above, however, such a large radiative annihilation cross-section would be too large to be compatible with naive expectations for a thermal relic, at least in our simple scenario and for standard cosmological and astrophysical assumptions. In any case, although the observation of a VIB or line-like feature from DM annihilation is a fascinating possibility, we caution that more data and a much more refined analysis, taking into account all systematics of the LAT, are required to reject or confirm this interpretation.
Under the above assumptions, the optimal target region is given by the set of pixels T o for which the overall signal-to-noise ratio R To , defined as is maximized. Finding the true T o requires in principle a scan over all ∼ 2 180 2 possible pixel combinations. Since this is unfeasible, we obtain an approximate T o by using the following simple algorithm: 1. We start with an empty set T and include only the one pixel with the largest individual signal-to-noise level R i as seed (this pixel typically lies at the galactic center).
2. For each pixel that is not already in T , we calculate how R T changes when this pixel is added; then all pixels that are found to improve R T are added to T at once.
3. For each pixel in T , we calculate how R T changes when this pixel is removed; then all pixels for which an increase of R T is found are removed from T at once. 4. We repeat steps 2 and 3 until the number of pixels in T remains constant.

5.
We remove remaining small regions in T that are not directly connected to the (always dominating) region at the galactic center.
The target region obtained in this way is a very good approximation to the optimal region, T o ≃ T . Note that the removal of remaining small regions in point 5 does not affect our results but merely cleans up the borders of the derived target region, and that the final regions are practically independent of the position of the seed pixel in point 1. The regions obtained in this way for the different adopted DM profiles are shown in Fig. 3, and will be used during our spectral analysis.

B Details on the statistical analysis
In order to study the sampling distribution of the test statistic TS (3.3) in absence of a signal, we performed a subsampling analysis of Fermi LAT data. To this end, we extracted the gamma-ray events measured in the hemisphere pointing towards the anti-galactic center, with longitudes |ℓ| > 90 • . Any signal from DM annihilation should be significantly suppressed in this direction. From these events we generate 30000 random sample spectra, with the Poisson expectation values in each energy bin given by µ j = f c j . Here, c j is the number of actually measured events in bin j, and f = 0.13 is adjusted such that the total number of events above 1 GeV in each random sample is ∼ 4 × 10 5 (in Reg1 and Reg2 the number of events are 5.8 × 10 5 and 2.7 × 10 5 , respectively). In the limit of large event numbers, this is equivalent to subsampling the energy distribution of the measured events with replacement.
In each of these sample spectra, we search for VIB features like discussed above and record the largest TS value found. The histogram of the maximal TS values that are obtained in this way is shown in Fig. 9. There, we also show the distribution that one obtains when selecting the maximum from 4 trials over a χ 2 k=2 distribution. The agreement is very good, and we used this distribution when calculating the look-elsewhere effect above. Histogram 4 trials over a χ 2 k =2 distribution Figure 9. Histogram of the maximal TS values obtained from a subsampling analysis of the Fermi LAT data in the hemisphere pointing towards the galactic anticenter. The blue line shows the theoretical distribution that we used to calculate the look-elsewhere effect.
In Fig. 10 we show the observed limits (black solid lines) in comparison with the limits that are expected at 68% (yellow) and 95% (green) CL. We derived these expected limits from 2000 mock data samples that were generated from the null model. In Reg1 to Reg3, the limits at m χ ≃ 150 GeV are significantly weaker than the expectation; this corresponds to the large TS values in the left panel of Fig. 5. On the other hand, the relatively strong limits at around m χ ≈ 100 GeV and m χ ≈ 250 GeV are a consequence of the 150 GeV excess, which influences the background fits. To illustrate this, we show by the dashed black lines the limits that we obtain when removing all data between 115 to 145 GeV (where the excess is most pronounced) from the fits; in this case the limits remain in the expected range.