ournal of C osmology and A stroparticle hysics J Model-independent indirect detection constraints on hidden sector dark matter

. If dark matter inhabits an expanded “hidden sector”, annihilations may proceed through sequential decays or multi-body ﬁnal states. We map out the potential signals and current constraints on such a framework in indirect searches, using a model-independent setup based on multi-step hierarchical cascade decays. While remaining agnostic to the details of the hidden sector model, our framework captures the generic broadening of the spectrum of secondary particles (photons, neutrinos, e + e − and ¯ pp ) relative to the case of direct annihilation to Standard Model particles. We explore how indirect constraints on dark matter annihilation limit the parameter space for such cascade/multi-particle decays. We investigate limits from the cosmic microwave background by Planck , the Fermi measurement of photons from the dwarf galaxies, and positron data from AMS-02. The presence of a hidden sector can change the constraints on the dark matter by up to an order of magnitude in either direction (although the eﬀect can be much smaller). We ﬁnd that generally the bound from the Fermi dwarfs is most constraining for annihilations to photon-rich ﬁnal states, while AMS-02 is most constraining for electron and muon ﬁnal states; however in certain instances the CMB bounds overtake both, due to their approximate independence on the details of the hidden sector cascade. We provide the full set of cascade spectra considered here as publicly available code with examples at http://web.mit.edu/lns/research/CascadeSpectra.html.


Introduction
Indirect searches provide one of the best ways to probe the nature of dark matter (DM) beyond gravitational interactions. Through the observation of gamma rays, cosmic rays, and the anisotropies of the Cosmic Microwave Background (CMB), we may find a hint of DM annihilations to Standard Model (SM) particles. Many models have been proposed in which DM annihilates directly to a pair of SM particles through, for example, a Higgs [1,2], gauge boson [3], axion [4], or neutrino [5]. Going beyond these simple models, we can consider scenarios in which DM is secluded in its own rich dark sector; such a setup is well motivated from top-down considerations (e.g. [6] and references therein). In such scenarios, the DM does not couple directly to SM particles (or such couplings are highly suppressed), but instead annihilates to unstable dark sector particles. These states may decay to SM particles or to other dark sector states, but eventually mediator particles that couple to the SM are produced. The mediators subsequently decay into SM particles, which in turn decay to stable and detectable photons, neutrinos, electrons, positrons, protons and/or antiprotons. We refer to this pattern as a "cascade annihilation" or simply "cascade", with a number of steps given  The DM, secluded in its own hidden sector, first annihilates to a pair of hidden sector particles. These φ n mediators subsequently decay to lighter particles in the hidden sector and finally to SM particles. Here we consider SM = {γ, e, µ, τ, b, H, W, g} and n = 0-6 step cascades (where n = 0 refers to the usual case of direct annihilations). Right: an equivalent diagram depicting the case where the DM annihilates through an off-shell heavy mediator; effectively decaying to an n-body state in the hidden sector which then decays to SM particles. by the number of distinct on-shell dark-sector states between the initial DM annihilation and the production of SM particles. We illustrate this setup schematically in figure 1.
Hidden Sector DM scenarios encompass a broad class of models. For instance models containing one light dark photon mediator [7][8][9], generically give rise to one-step cascades decays; DM annihilates to two dark photons which decay to SM particles. Multi-step cascades can occur naturally in hidden valley models [10,11]. In such models, production of the DM at terrestrial colliders and scattering in direct detection experiments can be generically suppressed by the small coupling between the dark and visible sectors. In contrast, indirect detection signals depend primarily on the annihilation rate of the DM to particles within the dark sector; the small coupling between the sectors only suppresses the decay rate of the mediators to SM particles, which does not affect indirect searches provided the decay rate is small on astrophysical timescales (as in e.g. [12]). Thus cascade annihilations scenarios are often invoked to explain anomalies and suggest new DM signals. For instance in [13][14][15][16][17][18][19][20][21] multistep cascades were used to explain the apparent excess GeV gamma-rays identified in the central Milky Way [22][23][24][25][26][27][28][29][30][31][32][33], while evading bounds from DM direct detection experiments. 1 In general the injection of photons and other high energy secondary particles produced is constrained by a number of indirect searches. In particular we focus on: • Measurements of the CMB by Planck [36] • Bounds set by Fermi from DM searches in the Dwarf Spheroidal Galaxies of the Milky Way [37] • Measurements of the e + flux by AMS-02 [38,39] Constraints from the above three experiments can be parametrized model-independently for the case of direct DM annihilations (see for instance [40]), by classifying annihilations to all possible two-body SM final states, DM + DM → SM + SM. For a given DM mass and final JCAP06(2016)024 state, the spectra of secondary particles, is fixed independently of the form of the DM interaction and spin. Therefore constraints on DM annihilation rates are usually quoted in terms of the parameters relevant to the direct annihilation scenario, and do not encompass DM models embedded in a hidden sector. 2 Given the broad space of Hidden Sector DM models, it is essential to provide model-independent methods that cover the majority of model space.
In the present work, we present DM mass dependent bounds on the DM cross section from the above three indirect detection experiments for DM annihilations via 0-6 step cascades to eight SM final states: γγ, e + e − , µ + µ − , τ + τ − , bb, gg, W + W − , and hh. We remain agnostic about the details of the hidden sector, thus making our statements robust and modelindependent. Limits from the Fermi dwarfs and AMS-02 generally provide the strongest robust constraints on channels that are rich in photons and those that are not, respectively (although at sufficiently high masses, limits from H.E.S.S. [42] and VERITAS [43] overtake those from Fermi ). While there may be arguably stronger bounds from the Galactic Center (e.g. [44]) or galaxy clusters (e.g. [45]), these limits depend strongly on the assumed DM density profile and/or the degree of substructure. We include the CMB limits because they are robust and almost independent of the spectrum of the annihilation products; thus we expect them to be nearly unaffected by the transition from 2-body to multi-body SM final states.
Our results are presented in figures [8][9][10][11], and our findings can be summarized as follows: • The Planck CMB bounds are robust and nearly model-independent varying by at most a factor of 1.5 over cascades with up to 6 steps for all final states.
• For photon-rich final states (all states considered except electrons and muons), we find the dwarf limits yield the most sensitive robust constraint, and can be weakened or strengthened by about an order of magnitude or more as compared to the direct annihilation case. For high (low) DM masses and small (large) step number the dwarf bounds can be overtaken by the robust CMB bounds as the most limiting constraints.
• For final states with few photons (electrons and muons), constraints from AMS-02 generally dominate the limits for low number of cascade steps. The limits can change by several orders of magnitude as compared to the direct case. As these weaken for higher DM masses and larger number of steps, CMB constraints become more important.
• Taking the above three points into account we find that for a fixed DM mass and final state, the presence of a hidden sector can change the overall cross section constraints by up to an order of magnitude in either direction (although the effect can be much smaller).
In addition to these constraints we also discuss how the bounds from multi-step cascades can be generalized to include the case of decays to n-body states in the dark sector. Finally as a supplement to this work we release code to generate the cascade spectrum.
In section 2 we review the procedure used in [13] to calculate the photon, electron and positron spectra from a multi-step cascade. Section 3 contains a description of the SM final state spectra used. Then in section 4 we describe how results for multi-body decays can be estimated from our cascade results. Our main results are presented in section 5-7 where we show the model-independent bounds extracted from the CMB, dwarfs, and AMS-02 respectively. We briefly discuss estimated constraints from antiproton data in section 8. Figure 2. The 0-step or direct photon (left), positron (center) or antiproton (right) spectrum for the various final states considered in this work. We have left out the γγ spectrum in the photon case and the electron spectrum in the positron case as both of these are δ-functions. Where applicable spectra are plotted with ǫ f = 0.3 or m φ = 20 GeV in the case of gluons.

JCAP06(2016)024
We discuss the interplay of the various experimental limits in section 9, present an example application of these results to the Galactic Center excess in section 10, and conclude in section 11. In the appendices we describe the contents of the publicly available code, as well as additional details and cross-checks.

Multi-step cascade annihilations
The multi-step cascade annihilation scenario is illustrated schematically in the left panel of figure 1. In this setup the DM pair annihilates into two scalar mediators (the case of non-scalar mediators was discussed in [13] where the conclusions proved to be relatively insensitive to choice of vector or scalar mediator 3 ) which subsequently decay through a (possibly) multi-step cascade in the dark sector, eventually producing a dark-sector state (with high multiplicity) that decays to the SM. Schematically we have: Here n is the number of steps as defined above. A variation on this picture occurs when any of the heavy hidden sector mediators goes off-shell and can therefore be integrated out yielding an effective vertex, now with a multibody decay in the hidden sector of the form φ n → mφ n−1 , with m > 2. This possibility is illustrated schematically in the right panel of figure 1 and for a 1-step cascade the analogue to eq. (2.1) would be: from which the extension to higher step cascades should be clear. Naively this framework seems like it could give quite different results to iterated 2-body decays. Yet in both cases, the main effect is to distribute the energy of the annihilation among a larger number of particles, thus increasing the multiplicity of the SM final state, lowering the average energy of the annihilation products, and broadening their spectrum. Consequently, limits on such scenarios can be broadly understood in terms of the n-step cascade results. This again highlights the point emphasized in [13] that the simple framework of n-step 2-body scalar 3 A thorough investigation of possible exceptions to this result is left to future work. cascades can describe a wide class of models and in this sense provide a relatively modelindependent framework. In eq. (2.1) and eq. (2.2) "(SM final state)" denotes the SM particles produced by a single decay of φ 1 , which in turn will (in general) subsequently decay to produce observable photons, neutrinos and charged stable particles. For example a SM final state may produce additional photons due to final state radiation (FSR) or the decay of neutral pions produced during hadronization. The mass ratio between φ 1 and the sum of the masses of the SM particles in this state, which we denote ǫ f (ǫ f ≡ ( m SM )/m 1 ), controls the level of FSR and hadronization, and so is a useful parameters for describing these decays; the details are discussed in [13]. When the SM particles are massless, the relevant parameter is instead just the mass of the φ 1 , which we denote interchangeably as m 1 or m φ .
The spectrum of particles in an intermediate step of a cascade may be obtained using the method discussed in [13], which we briefly review in this section. Consider the "ith step" decay φ i+1 → φ i φ i . In the rest frame of φ i+1 we will denote the spectrum of the subsequent photons, electrons or positrons as dN/dx i , where and E i is the energy of the photon, electron or positron in the φ i+1 rest frame. We define ǫ i = 2m i /m i+1 , and will (by default) assume a large mass hierarchy between cascades steps such that ǫ i ≪ 1. Assume that the spectrum in the rest frame of the φ i particle is known and denoted by dN/dx i−1 . In the limit of a large mass hierarchy the decay of φ i+1 produces two highly relativistic φ i particles, each (in the rest frame of the φ i+1 ) carrying energy equal to m i+1 /2 = m i /ǫ i . The photon, electron, or positron spectrum per annihilation in the rest frame of the φ i+1 is then given by a Lorentz boost, and takes the simple form In this way, we can begin with a direct spectrum of dN/dx 0 from φ 1 → SM final state -the details of which are described in the next section -and generate a cascade spectrum inductively. By repeated application of this formula we can see that the presence of each additional step in a cascade acts to broaden and soften the spectrum, and shift the peak to lower masses. Importantly the shapes of these cascade spectra are very simple, being characterized by just three pieces of information: the number of steps n, the SM final state (often denoted f ), and the value of ǫ f . Such cascades are independent of the details of each of the intermediate steps, within the large-hierarchy (ǫ i ≪ 1) approximation, and as such are independent of the various ǫ i . 4 As pointed out in [13], although the large-hierarchy approximation seems to discard information, the more general case can be recovered quite easily. To see this, consider the opposite limit where ǫ i → 1, so that 2m i ≈ m i+1 . In this case, the rest frames of the φ i+1 and φ i are the same, so no boost needs to be applied. As such, in this "degenerate limit", the final spectrum of annihilation products is the same as that for a hierarchical cascade with one fewer step, with half the initial DM mass and half the annihilation cross-section. The intermediate regime, where neither ǫ i nor 1 − ǫ i are particularly small, smoothly interpolates between these two cases. Thus by studying the parameter space of (m χ , σv , no. of steps) in the hierarchical limit, it is possible to quickly estimate results for a general cascade.
Again this framework is more general than it might initially appear. For example, simple extensions where a φ i decays to two φ i−1 with different masses will not change our results

JCAP06(2016)024
in the large-hierarchy limit, as those results are independent of the intermediate masses.
Additionally, as pointed out in [13], for larger n our cascade scenarios can approximate models with hadronization in the dark sector (see e.g. [46,47]), and additionally as we will show in section 4, multi-body decays can also be approximately captured within this framework.
Note that the cascade scenario must be physically self-consistent: the mass hierarchy between the DM mass and the SM particles in the final state must be sufficiently large to accommodate the specified number of steps. In other words, there is a hard upper limit on the number of steps allowed, for a given DM mass and final state. In detail, for an n-step cascade ending in a final state consisting of two particles each with mass m f , we defined ǫ f = 2m f /m 1 , ǫ 1 = 2m 1 /m 2 , ǫ 2 = 2m 2 /m 3 and so on until ǫ n = m n /m χ . Combining these, the DM mass is given in terms of m f and the ǫ factors by: If the ǫ i factors are allowed to float, we can still say that 0 < ǫ i ≤ 1 in all cases (since each decaying particle must have enough mass to provide the rest masses of the decay products), setting a strict lower bound on the DM mass of: Where this limit is not satisfied, the spectra should not be thought of as potentially representing a physical dark-sector scenario, but only as a parameterization for general spectral broadening. For the massless final states considered here (photons and gluons) m f = 0, but we can still derive a condition from the value of m φ , specifically: 3 Direct spectra Using the formalism outlined in the previous section, from a given "direct" spectrum we can straightforwardly generate an n-step cascade spectrum, to compare with various indirect searches. We outline the different SM final states considered for the direct (0-step) spectra in this section. To obtain limits using bounds from the dwarfs, CMB and AMS-02 we need the spectrum of photons, electrons and positrons, and so we determine the spectrum for these particles arising from the boosted decays of the following eight SM states: γγ, e + e − , µ + µ − , τ + τ − ,bb, W + W − , hh, and gg. 5 We choose these states as a representative sample of possible spectra. For example decays of light quarks generally give signals similar to those of b-quarks and the ZZ final state is similar to W + W − . As discussed in the previous section, many of the cascade spectra depend on the parameter ǫ f = m SM /m 1 = 2m f /m 1 (the final equality holds for all the processes we consider here). In the context of generating the direct (0-step) spectrum, we can imagine two analogous processes: either the direct annihilation χχ → SM final state, in which case ǫ f = m f /m χ , or the final step in a cascade annihilation, φ 1 → SM final state, so that ǫ f = 2m f /m 1 as stated. If the (SM final state) is a photon or a gluon, then clearly ǫ f is no longer a useful parameter; instead m φ = m 1 (equivalent to 2m χ in the case of direct annihilation) plays this role. For JCAP06(2016)024 many spectra no such parameter is needed. For example the γγ photon spectrum, as well as the positron spectra from e + e − or µ + µ − final states, are independent of any such parameter, since they are either just δ-functions or arise from decay rather than FSR or hadronization.
In all but five cases, we use the results of the PPPC4DMID package [48] to produce the direct spectra (hereafter referred to simply as PPPC). The exceptions to this are: • the γγ photon and e + e − electron or positron spectra, which are just δ-functions, to a good approximation (we neglect the effect of FSR on the e + e − spectra in the case of annihilation/decay to e + e − ), • the spectra of photons produced in conjunction with the e + e − and µ + µ − final states, for which we use the analytic results of [13,49] and [13,49,50] respectively, • the spectrum of electrons or positrons from muon decay, where we use the analytic Michel spectrum [49,51].
Finally we briefly comment on the ǫ f or m φ dependence of the various direct spectra as it is often useful in interpreting results, noting that [13] has a more detailed discussion of several cases for photon spectra. For photons produced from e + e − and µ + µ − final states, the spectra arise entirely from FSR and so are strongly dependent on ǫ f , increasing in flux and becoming more sharply peaked near the maximum energy as ǫ f decreases. Similarly the photon spectrum produced from the W -boson final state, in addition to a broad continuum peaked at low x, acquires a sharp spike at high x due to FSR when ǫ f becomes small. The photon spectrum from the b-quarks final state broadens and moves its peak to smaller x as ǫ f decreases; the gluon spectrum behaves similarly as m φ increases. Finally the photon spectra from τ + τ − andhh final states are largely independent of ǫ f .
The positron spectra produced from the Higgs and tau final states again show no real variation with ǫ f . For positrons the spectrum from the W + W − final state is also quite independent of ǫ f , whilst the b-quark and gluon spectra behave much as they did in the photon case. Lastly, for antiprotons, once more the spectra from Higgs and W -boson final states are independent of ǫ f , whilst now for decreasing ǫ f (increasing m φ ) the b-quark (gluon) spectrum increases in height without substantially changing the position of its peak.

Multi-body cascades
So far we have focused on cascades comprised of 2-body scalar decays. In this section we discuss the extension of this framework to the case of n-body cascades, schematically illustrated on the right of figure 1. As we will see, in the large hierarchies regime the n-body decays can be understood in terms of our existing 2-body results, again emphasizing the modelindependence of our results. The explicit calculations and examples to help build intuition are provided in appendix A.
As explained in the introduction, a multi-body decay can arise if there is a heavy mediator in the cascade that has been integrated out. This can happen anywhere in a cascade, but here we restrict to a 1-step cascade of the form χχ → n × φ → 2n × (SM final state) (cf. eq. (2.1)). From here the extension to higher step cascades is intuitively clear, and in practice can be calculated using eq. (2.3). As shown in the appendix, an analogue of this equation can be derived for the multi-body case: Figure 3. The 1-step spectrum for an n-body cascade from a direct annihilation to bb with ǫ b = 0.1 is shown as the dashed gray curves for n = 2-10, where lighter curves correspond to larger n. In purple, green and orange we show a 2-body 1-step, 3-step and 5-step cascade spectrum respectively, for the same direct spectrum. These three curves outline the n-body results and show that the result of 1-step multi-body spectra should be encapsulated in the multi-step 2-body results.
where again dN/dx 0 represents the direct spectrum. Intuitively, the dx 0 integral accounts for the boosting of the decay products, just as in eq. (2.3), whilst the ξ integral samples from the n-body phase space to give the correct degree of boosting. At first glance it appears that this formula could produce marked differences to our standard cascade framework, but as we show in figure 3, this is not the case. There we show the 1-step spectrum for an n-body cascade ending in annihilation into the SM state bb with ǫ b = 0.1, for n = 2-10. Overlaid is the m-step 2-body cascade for m = 1,3,5. Of course when n = 2 and m = 1 these two are the same by definition. But increasing n and increasing m perturb the spectra in quite similar ways (albeit to different degrees, as expected since the multiplicities of final-state particles are not equal for m = n with n > 2), and so we expect the observational signatures of multi-body decays to lie within the space mapped out by the simple cascade annihilations. An example of the constraints on multi-body decays, and how they lie within the band of cascade constraints, is given in figure 14, in appendix A.

CMB bounds from Planck
DM annihilation during the cosmic dark ages can inject ionizing particles into the universe, modify the ionization history of the hydrogen and helium gas, and consequently perturb the anisotropies of the CMB. Sensitive measurements of the CMB by Planck [36] (and previously WMAP and other experiments) can place quite model-independent limits on such energy injections, which when applied to DM annihilation are competitive with other indirect searches.
The figure of merit for CMB limits on DM annihilation is the parameter p ann = f eff σv /m χ , where f eff is an efficiency factor that depends on the spectrum of injected electrons and photons, and the other factors describe the total power injected by DM annihilation per unit time. In principle, different DM models could give rise to different patterns of anisotropies in the CMB -but for WIMP models of DM that annihilate through s-wave processes, it has been shown [52] that the impact on the CMB is identical at the sub-percent -8 -

JCAP06(2016)024
level up to an overall rescaling by p ann (related studies [53][54][55] independently found that the signal was largely controlled by a single parameter). In [56]- [57] this result was generalized to include any class of DM models for which σv can be treated as a constant during the cosmic dark ages, which is generically true for the models considered in the present work. We use the results of [56] to compute f eff as a function of DM mass and annihilation channel. In particular, we compute positron and photon spectra for the case of direct annihilations using PPPC, and then convolve to find the resulting spectrum for an n-step cascade as discussed above. The spectrum of electrons is equal to that of positrons by the assumption of charge symmetry. We then integrate the resulting spectra over the f eff (E) curves provided in [56] to obtain the weighted f eff (m χ ) for n = 0-6 step cascades for DM annihilations to various final states: We neglect the contribution from protons and antiprotons, as for all the channels we consider, the fraction of power proceeding into these channels is rather small, and consequently including them should only slightly increase f eff [58]. The constraints we present are therefore somewhat conservative (they could be strengthened slightly by a careful treatment of protons and antiprotons). As discussed in [57], we use the best-estimate curves suited for the "3 keV" baseline prescription, which are most appropriate for applying constraints derived by Planck.
The bound set on the annihilation parameter, p ann , from Planck temperature and polarization data is taken to be [36]: In figure 8 we present our results for the bound on DM cross-section as a function of m χ for various numbers of cascade steps and SM final states. We note that the number of steps does not affect the total power deposited by DM annihilation per unit time (at least in the simple scenario where all that power eventually goes into SM particles). Each additional step reduces the average energy of the final-state photons/positrons/electrons by a factor of 2, but simultaneously increases their multiplicity by a factor of 2. Thus the only possible impact on the constraints comes from the energy dependence of f eff , combined with the softening and broadening of the spectrum.
In accordance with our expectations, we find that the effect of the spectral broadening and softening is rather mild, typically changing the constraints by no more then 0.1-0.15 decades (corresponding to a factor of ∼ 1.5). There is no general trend, in that constraints on these high-multiplicity final states may be either weaker or stronger than those pertaining to direct annihilation; this arises from the fact that f eff is not a monotonic function of energy, so lowering the average energy of the injected particles may either increase or decrease the deposition efficiency. In general, f eff and hence the upper bound on the ratio σv /m χ varies less as a function of mass for higher-multiplicity final states (as expected, from the broader resulting spectrum), but this effect is very small. The choice of ǫ parameters, again, does not perturb the constraints outside this ∼ 0.15-decade band. We refer the reader to the appendix B for additional details regarding the behavior of f eff . The dwarf spheroidal galaxies of the Milky Way are expected to produce some of the brightest signals of DM annihilation on the sky. Whilst less intense than the emission expected from the galactic center, the dwarfs have the distinct advantage of an enormous reduction in the expected astrophysical background. These features make them ideal candidates for analysis with the data available from the Fermi Gamma-Ray Space Telescope. Indeed the Fermi Collaboration has set stringent limits on the DM annihilation cross-section using the dwarfs [37], and together with the DES Collaboration have used 8 newly discovered dwarf satellites [59,60] to set independent limits [61]. We note in passing that several groups have pointed out an apparent gamma-ray excess in the direction of one of the new dwarfs, Reticulum II [61][62][63], albeit with considerable variation as to its significance (with estimates ranging from ∼ 3σ to completely insignificant). We will not discuss this tentative excess here, other than to note as it appears roughly consistent with the emission coming from the GCE, the implications for dark sector cascades will be analogous to those discussed in [13].
Here we focus on understanding how the presence of cascade annihilations can modify the limits obtained from these dwarf galaxies. In order to do this we use the publicly released bin-by-bin likelihoods provided for each of the dwarfs considered in [37]. 6 This analysis made use of 6 years of Pass 8 data and found no evidence for an excess over the expected background. Note the Fermi collaboration produced an earlier analysis of the same dwarfs using 4 years of Pass 7 data in [64]. In appendix C we show that the results are similar between the two, but that the limits set using the newer analysis are usually about half an order of magnitude stronger.
For a given dwarf Fermi provides the likelihood curves as a function of the integrated energy flux in each of the energy bins considered in their analysis, covering the energy range from 500 MeV to 500 GeV. Thus to obtain the likelihood curves for our cascade models we need to firstly determine the integrated energy flux per bin. This will be a function of the DM mass m χ , annihilation cross-section σv , and shape of the cascade spectrum dN/dx -which itself depends on the number of cascade steps, the identity of the final state particle and possibly either ǫ f or m φ . For an energy bin running from E min to E max , the energy flux in GeV/cm 2 /s is: where J i is the J-factor appropriate for the individual dwarf i. Treating the energy bins as independent, we can simply multiply the likelihoods for the various bins to obtain the full likelihood for a given dwarf i: L i (µ|D i ), which is a function of both the model parameters µ and the data D i . At a given mass and for a given channel (final state and number of cascade steps), µ just describes the annihilation cross-section σv . There is, however, one additional source of error that should be accounted for: the uncertainty in the J-factor. Following [37] we incorporate this as a nuisance parameter on the global likelihood, modifying the likelihood where for log 10 (J i ) and σ i we use the values provided in [37] for a Navarro-Frenk-White profile [65]. This approach allows us to account for the J-factor uncertainties using the profile likelihood method [66]. We obtain the full likelihood function by multiplying the likelihoods for each of the 15 dwarfs together. Using this likelihood function, for a given DM mass and cascade spectrum we can then determine the 95% confidence bound on the annihilation cross-section. We follow this procedure for cascade annihilations with 0-6 steps, for final state electrons, muons, taus, b-quarks, W -bosons, Higgses, photons and gluons, considering two different values of ǫ f or m φ where appropriate.
Results are shown in figure 9. For the final states considered in [37], our direct/0-step results are in agreement. Recall that there is a physical limitation on realizing a given cascade scenario set by m χ ≥ 2 n m f /ǫ f , as mentioned in section 2. The constraints corresponding to scenarios that satisfy this condition are indicated by darker lines, but we also show the limits for cases that do not satisfy this condition (and so cannot be physically realized as a cascade annihilation of the type we have considered), to demonstrate the effect of spectral broadening.
Before discussing results for each final state independently, there are a few generic features worth pointing out. Recall that higher-step cascades have a spectrum peaked at lower x = E γ /m χ . Thus in order to produce emission at an equivalent energy, higher-step cascades require a larger DM mass, which in turn requires a larger cross-section to inject the same amount of power (as the DM number density scales inversely with the mass). Equivalently, at a fixed mass and cross-section, larger numbers of cascade steps will tend to produce a larger number of lower-energy photons; at low masses, some of these photons may lie outside the energy range of the Fermi analysis, and the astrophysical backgrounds will also generally be larger at low energies. These factors tend to weaken the constraints, and indeed we see a systematic trend for weaker bounds with increasing n for low-mass DM, for all channels.
Nevertheless this conclusion is not inevitable. Specific energy bins may allow stronger constraints than neighboring bins, purely due to statistical accidents; adding cascade steps smooths out such effects. The total number of emitted photons is increased with larger n (albeit while preserving the total injected power).
Most generically, if the DM mass is large, much of the spectrum may be above the 500 GeV cutoff of this analysis in the case of direct annihilation. In this case, adding cascade steps can strengthen the constraints by moving the photons into the range of sensitivity for the search. This effect is most pronounced, and occurs at the lowest DM masses, for final states with spectra peaked at large x (electrons, muons, taus and photons): for softer directannihilation spectra, even at the heaviest masses tested, the peak of the spectrum does not move past 500 GeV. Inclusion of higher-energy data, e.g. from studies of the dwarf galaxies with VERITAS [43], would potentially strengthen the constraints at high DM masses, but for this reason we expect the improvement to be smaller for higher-step cascades.
Thus in general we see a weakening in the cascade constraints relative to the directannihilation case at low DM masses, and a strengthening at high DM masses, with the crossover point and the width of the band varying based on the SM final state. For some final states, the cascade constraints can be weaker or stronger than those for the direct- annihilation case by more than an order of magnitude. Let us now discuss the detailed results for each SM final state (shown in figure 9) separately: Electrons: the generic behaviors discussed above are clearly demonstrated in these results. There is also a striking difference between the results for direct and cascade annihilations. The photon spectrum in the direct case originates from FSR and is very sharply peaked (especially for small ǫ f ); even a single cascade step will smooth out the spectrum and considerably change its shape. Further, the bounds are strongly dependent on the value of ǫ f , as the FSR photon spectrum diverges as ǫ f → 0. As such, for smaller ǫ f we expect stronger limits, and this is exactly what we observe. Nonetheless note that the position of the peak of the spectrum in x is not strongly dependent on ǫ f , so we should expect the crossover behavior between different spectra mentioned above should happen at a similar location for different ǫ f values and this is exactly what we observe. Finally note that the bumps in the direct spectrum are a result of the sharply peaked 0-step spectrum moving between energy bins. The width of these bumps is exactly the width of the energy bins in the data. As we move to cascade scenarios, the spectrum is smoothed out and the majority of the emission is no longer in a single bin, meaning these bumps vanish.
γγ: the most noticeable feature here is the jagged direct spectrum. As the direct spectrum of γγ is just a δ-function at the mass considered, these jumps are an extreme realization of the issue mentioned for the 0-step electron limits: we get a jump as the emission moves from one of the energy bins considered to the next. Of course physically the Fermi instrument has a finite energy resolution, which will act to smooth out such a sharp feature. To approximate this we smooth the 0-step by a Gaussian with a width set to 10% of the energy value, yet this ultimately had little impact on the extracted limit. Note also that once the spike moves beyond 500 GeV, which occurs at roughly log 10 m χ = 2.7, the Fermi data can no longer constrain this scenario so the limit completely drops off.
Muons: the photon spectrum for the muon final state is very similar to that for the electron final state, except that it is slightly less dependent on ǫ f . The results here are accordingly very similar to those for the electron final state, except that the variations with ǫ f are less pronounced.
Taus: the fact that the tau spectrum is only weakly dependent on ǫ f is clearly visible; otherwise only the generic behavior is apparent. b-quarks: there is a modest dependence on ǫ f , which does not change the qualitative results. The crossover where the direct constraints become weaker than the cascade constraints occurs at a DM mass around 100 GeV. Due to the kinematic bounds, over the physically allowed region the variation in the band width is fairly modest varying by at most 0.4 decades.
Gluons: the gluon spectrum behaves very similarly to the b-quark spectrum, if we swap decreasing ǫ f for increasing m φ . As such the results are similar to those for b-quarks.
W -bosons: firstly note that the kinematic edge in these results appears from the threshold requirement to have enough energy to create on-shell W 's. Other than this we see that the limits are somewhat stronger for smaller values of ǫ f , which is because the W spectrum includes a small FSR component which is larger for smaller ǫ f . The width of the band of possible results is at most 0.7 decades. Again we also see a crossover where the direct constraints become weaker than the cascade constraints, here at roughly 500 GeV. Higgs: as with the W -bosons, our results again have a kinematic edge. Furthermore, like final state taus, the Higgs spectrum is only weakly dependent on ǫ f and thus so are the results. As for the W case, the width again has a maximum around 0.7 decades, whilst this time the direct crossover first happens at about a TeV.

Positron bounds from AMS-02
AMS-02 has recently released a precise measurement of cosmic ray electrons and positrons in the energy range of ∼ 1 GeV to ∼ 500 GeV [38,39]. The measured positron ratio exceeds the prediction of the standard cosmic ray propagation models at energies larger than ∼ 10 GeV. There are many possible explanations for this rise in the positron ratio, including DM physics (although the annihilation scenario seems challenged by a range of other null results, e.g. [36]), nearby pulsars [67,68] or supernovae [69].
The presence of an apparent large positron excess of unknown origin makes it challenging to set stringent limits on general DM annihilation scenarios. The situation is further complicated by the effects of solar modulation at energies below ∼ 10 GeV [70][71][72], which modifies the cosmic ray flux in a charge-dependent manner and adds significant astrophysical uncertainties. However, the data do indicate that both the positron ratio (flux of e + divided by the flux of e + + e − ) and the fluxes of cosmic ray electrons and positrons are fairly smooth; there is no clear structure in the spectrum within the energy resolution of AMS-02. Accordingly, it is possible to set quite strong constraints on DM models that predict a sharp spectral feature in the positron spectrum (e.g. [73][74][75]).
As discussed in the previous sections, DM annihilation through multi-step cascades usually gives rise to a softer and broader spectrum than direct annihilation to the SM states, generally leading to weaker bounds from AMS-02. In this section, we study this effect quantitatively. We note that our goal here is to study the impact of these spectral changes, not to explore possible explanations for the rise in the positron fraction or systematic uncertainties in the modeling of the background or signal.
To set bounds on annihilating DM, we first need to parametrize or model the backgrounds. Here the backgrounds that we refer to are the astrophysical cosmic ray flux, plus some new smooth ingredient to account for the observed rise in the positron flux. Since we do not know the source of the new ingredient, polynomial functions of degree 6 are introduced to fit the AMS-02 positron flux data (the 6 degrees are employed to obtain a good χ 2 fit to the data). To derive the limits, we float the 6 parameters from the polynomial functions within 30% of the best fit values from the fit without DM, together with the DM annihilation cross-section. We check that increasing the range of allowed values for the background parameters does not weaken the constraints.
We derive limits from only the positron flux, as both the positron and electron backgrounds are required to float in the fit to the positron ratio. Such an analysis would require many additional free parameters, and is beyond the scope of the current work. As a crosscheck we attempted a simplified fit to the positron ratio data (using AMS-02 measurements of the total e + + e − spectrum) and found constraints of comparable strength to those we present here.
The positron flux from DM annihilation is obtained by propagating the injected positron spectrum using the public code DRAGON [76,77]. There are substantial systematic uncertainties in the propagation of electrons and positrons in the galaxy, affecting diffusion, energy loss, convection, and solar modulation. In particular, accounting for uncertainties in Once electrons and positrons are injected into the halo, they will diffuse and lose energy. Their number density N i evolves according to the following diffusion equation, where D is the spatial diffusion coefficient, depending on the spatial position and energy. It is parametrized by the following form where we assume the diffusion zone is axisymmetric, and use the cylindrical coordinate system (r, z). Most of the electrons and positrons are trapped in the diffusion zone with thickness 2z t .
Here ρ = p/(Ze) is the rigidity of the charged particle with Z = 1 for electron and positron. D 0 normalizes the diffusion at the rigidity ρ 0 = 4 GV. In eq. (7.1), v c is the velocity of the convection winds;ṗ accounts for the energy loss; Q i is the source of the cosmic ray, where DM is one kind of the source; D pp accounts for the diffusion in the momentum space; the last two terms in eq. (7.1) parameterize how the nuclei inelastic scattering with the gas to affects the number density of the cosmic rays. Although there are many parameters in the diffusion equation, we do not simulate the backgrounds (instead modeling them with a polynomial function), which decreases the systematic uncertainties of the limit substantially.
We use a specific model to propagate the electrons and positrons injected by DM annihilation [78]. In this model, D 0 = 2.7 × 10 28 cm 2 /s, z t = 4 kpc, δ = 0.6 and we take the local density of DM to be 0.4 GeV/cm 3 and the density profile to be the Navarro-Frenk-White profile. We set the convection and diffusion in momentum space equal to zero, since they do not change the spectrum significantly in the energy range of interest [79]. There are other propagation models with different diffusion terms or diffusion zone heights that can be employed here. However, since the energy loss effect is dominant for the propagation of high energy leptons, we choose only one propagation model to derive the limits. While there may be remaining systematic effects due to the choice of propagation model, we reiterate that the purpose of this analysis is not to explore all the uncertainties in these constraints.
Cosmic-ray propagation is affected by the magnetic field, which determines both how the cosmic rays diffuse and their energy losses due to synchrotron radiation. The magnetic field is modeled by two components, one regular and one turbulent (irregular) [80,81]. [82] gives the constraints on these components. To be conservative, here we set the value of the magnetic field at the Sun to B ⊙ ∼ 8.9 µG. With this magnetic field, the local radiation field and magnetic field energy density is 3.1 eV/cm 3 , which is close to (but somewhat higher than) the 2.6 eV/cm 3 value used for conservative constraints in [73]. For this reason, the constraint we obtain for the direct annihilation is slightly weaker than even the conservative case studied in [73], as the energy loss rate for the positrons is higher. The main effect of changing the local energy density is to rescale all the constraint curves, with lesser effects on the variation of the constraint with DM mass and number of cascade steps.  For cosmic rays with energy smaller than 10 GeV, although there are many other parameters in the propagation model, we only consider the uncertainty from the solar modulation, which is modeled by the modulation potential. The modulation potential φ in the range of (0, 1) GeV is fixed by minimizing the χ 2 to fit the AMS-02 data.
In summary, we derive the limits on DM annihilation by using AMS-02 positron flux starting from 1GeV. The background is parametrized by a polynomial function of 6 degree, and to derive the bounds we let the 6 parameters float within 30% of their best fit values. The diffusion model is employed here to propagate the DM positron flux, and the solar modulation potential is allowed to float in the range (0, 1) GeV when minimizing the likelihood function. The limits are summarized in figure 10.
In general, similar to the dwarf galaxies, the constraints on cascade models can be substantially weaker than those on the direct-annihilation case for low DM masses (below ∼ 100 GeV), by up to several orders of magnitude depending on the channel. This weakening likely arises from a combination of (a) positrons falling below the minimum energy of the search, and (b) broadening of the spectrum making it easier for the background model to compensate for a DM component. The effect can be up to two orders of magnitude in most channels at sufficiently low masses (the exceptions are the W , Higgs and b-quark final states where low DM masses are kinematically forbidden). At high masses, the bands of possible constraints are narrower, of order half an order of magnitude or less; for the electron, muon, tau and gamma final states the direct-annihilation constraints are systematically weaker than those for cascade scenarios. This is likely due to the cascade scenarios producing greater -15 -

JCAP06(2016)024
numbers of positrons in the energy range of the search, but may also be due to the hardening of the positron flux at high energies mimicking a hard signal from DM annihilation.

Antiproton ratio bounds
Measurements of cosmic ray antiprotons may also provide stringent constraints on DM annihilation, especially to hadronic final states (e.g. [83][84][85]). However, as with the positron data, uncertainties in cosmic ray propagation and the modeling of the background cosmic ray population can substantially affect sensitivity. Here we present a brief discussion of the dependence of antiproton bounds on the number of cascade steps; we caution that this analysis does not include a comprehensive study of uncertainties in the background and propagation models, although we do study two alternative propagation models. Since it is difficult for antiproton limits to robustly improve on constraints from Fermi observations of dwarf galaxies (unlike the positron bounds, which have sensitivity to channels that produce very few photons), we do not include these limits in our final combined constraint plots. In this section we employ data from BESS [86,87], CAPRICE [88], and PAMELA [89]. 7 The propagation of antiprotons and protons follows the transport equation, as given in eq. (7.1). The parameters in the propagation models are derived by reproducing the B/C data, following the approach in [71,78,[90][91][92]. In this section, we will use the KRA and THN models [91] as illustrative examples of possible propagation scenarios; the THN model yields more conservative constraints on DM annihilation, since it possesses a thin diffusion zone (z t = 0.5) and small diffusion coefficient (D 0 = 0.32 × 10 28 cm 2 /s). We assume that the DM density is well described by the Navarro-Frenk-White profile, and that the local DM density is 0.4 GeV/cm 3 .
Instead of parametrizing the astrophysical background by polynomial functions, the astrophysical background of proton and antiproton are simulated by DRAGON. Besides the propagation parameters, we must also determine the solar modulation parameter φ. Since the different experiments whose data we use operated at different times, each one must be allowed to have a different value of this parameter. We determine φ for each experiment by minimizing the χ 2 considering only the ratiop/p; we then re-fit the data including a DM component but holding φ fixed.
As in section 7, we employ DRAGON [76,77] to propagate the flux from DM annihilation. Following [91], we derive our bounds from the antiproton-to-proton ratio data, to reduce systematic uncertainties and allow combination of results from different experiments.
The constraints on the DM annihilation cross section, for final states comprised of b quarks, are summarized in figure 5. For this figure we set ǫ f = 0.3; we also include the case of direct annihilation. We see that the choice of propagation model changes the constraints by up to an order of magnitude, with the constraints from the THN model being weaker.
We find that increasing the cascade step number strengthens these antiproton bounds for DM masses larger than ∼ 20GeV. The reason is that thep/p data are considerably more precise (and abundant) at energies smaller than 10 GeV; increasing the number of steps, for a given DM mass, shifts the antiproton spectrum down in energy, and into this better-constrained region. 7 Preliminary antiproton data have been presented by AMS-02, but we do not employ these results; we reiterate that the purpose of this section is not to provide a robust new bound on DM annihilation from antiprotons (as this requires an in-depth study of systematic uncertainties that is beyond the scope of this work), but to investigate broadly how these bounds may be expected to vary for multi-step cascades.  It is also worth noting that in this case, the direct annihilation bound does not behave like a limiting case of the bounds on an n-step cascade. This is because by holding ǫ f fixed for the n-step cascade, we are fixing the center-of-mass energy for the b quark production (i.e. the decay of the last mediator in the chain); in contrast, for the direct annihilation this quantity is determined by the DM mass. The antiproton production is sensitive to the details of the hadronization process and requires a relatively large center-of-mass energy; consequently the antiproton spectrum is quite dependent on the center-of-mass energy for the final decay/annihilation, in contrast to the electron and photon spectra. This suggests that our model-independent cascade formalism is likely to be a worse approximation for antiproton final states, with generically stronger dependence on ǫ f , compared to the other channels we have considered.
Nonetheless, in the case we have studied, even the conservative THN model admits constraints on high-step cascades that are stronger than those from the Fermi dwarf observations for relatively light DM (e.g. for DM masses below ∼ 100 GeV, for n ≥ 4). This suggests that improved antiproton data and better characterization of the background and cosmic ray propagation could be particularly important in constraining cascade-type models with hadronic final states.

General discussion
We summarize our main results in figure 11, where we overlay the combined constraints from the three experiments as a function of DM mass for an n = 0-6 step hidden sector cascade. Furthermore in figure 4 we show results just for the direct, or 0-step, annihilation, in order to highlight the interplay between the experiments. As discussed above the CMB constraint is fairly insensitive to the SM final state and number of cascade steps. The AMS-02 bounds, which are most constraining for direct annihilation for electron and muon final states, weaken rapidly at low masses as the positron spectrum broadens with increasing cascade steps, but for masses above a few hundred GeV are also fairly robust. The dwarfs are generally most constraining for final states with a high multiplicity of photons. However, for lepton-rich or photon-rich final states respectively, the CMB bounds can become more constraining than the -17 -

JCAP06(2016)024
AMS-02 or Fermi limits for large numbers of cascade steps at low masses, or small numbers of cascade steps at high masses. We summarize the results the various SM final states below.
Electrons: the spectrum of positron and photon spectrum is very sharply peaked in the case of direct DM annihilations to e + e − . Thus AMS-02 places the most constraining bound for n = 0-4 step cascades at low masses m χ 400GeV. As the number of steps increases, n > 3, the spectrum smooths and broadens thereby weakening the AMS-02 bound so that the CMB bound becomes the most constraining. The CMB bounds are generically stronger at high DM masses, above a few hundred GeV. The dwarf limits are, in all cases, 1-2 orders of magnitude less constraining than the AMS-02 and CMB bounds. This is unsurprising given the dwarfs are only sensitive to the photon spectrum from the final state electrons, which represents only a small fraction of the available power per annihilation.
γγ: the strongest constraints almost always arise from the Fermi dwarfs, although at high DM masses and for small numbers of steps, the CMB bounds may be more stringent. However, in this case VERITAS or H.E.S.S dwarf searches may actually provide a stronger limit. For AMS-02 the positron spectrum is similar in shape to that of the electron channel; the photon generates a hard electron spectrum via Drell-Yan. Nonetheless this process is suppressed by a factor of α e as well as phase space. Combining these, approximately two order magnitude suppression relative to electron case would be expected and is in fact observed.
Muons: recall that the spectrum of positrons and photons from DM annihilations to muons is similar to the corresponding spectra for the electron final state, except the photon spectrum is less dependent on ǫ f and the positron spectrum is somewhat broader. For 0-2 cascade steps, the most stringent constraints are from AMS-02 at low masses, below a few hundred GeV. At higher step numbers (for all masses) and higher masses (for all cascade scenarios), the CMB limit becomes more restrictive.
Taus: the tau final state is richer in photons than the other leptonic final states, and yields smoother and broader photon and positron spectra even in case of direct annihilation. Thus the bound from the dwarfs is more sensitive and constraining than AMS-02, and generally also stronger than the CMB limits. The exceptions are at low mass and large number of steps, or inversely high mass and a small number of steps, as in both cases the CMB bounds dominate the constraint.
b-quarks: the direct spectrum for DM annihilations to bb is much softer then the previously discussed channels. So the Fermi dwarf limits almost always provide the strongest constraint; for low masses and n = 3-6 steps there is a region of parameter space where the CMB bounds appear to be more stringent, however this region is kinematically disallowed.
Gluons: as previously discussed the gluon spectrum behaves very similarly to the b-quark spectrum, if we swap the decreasing ǫ f for increasing m φ . As such the results are similar to those for bb.
W -bosons: for annihilation to W final states the bounds are quite robust, with the dwarfs always setting the strongest limits.

JCAP06(2016)024
Higgs: annihilations to final state Higgses is similar to the W case; the results are almost identical aside from the difference in the kinematic edge between the H and W mass.
Finally let us say a few words about how these results are likely to change in the near future as additional experimental data becomes available. For the CMB limits, the shape of the limits is dictated by processes occurring in the early universe -additional data from Planck will strengthen the limits but not change their qualitative behavior. The situation is similar for the dwarf limits. In the absence of a signal the shape of our results are being determined by three factors: the DM spectrum, the spectrum of the backgrounds, and the response of the Fermi detector. These will not change with additional data. 8 Again additional data and especially the addition of new dwarfs will just strengthen the limits (unless, of course, there is a detection). Given AMS-02 is still relatively early in its mission cycle, it is possible the relative sensitivity at different energies could change as more data is collected, which could alter the behavior of the constraints with respect to the number of cascade steps; however, the main effect is likely to be just to strengthen the limits. The addition of AMS-02 antiproton data to our estimated antiproton bounds could likewise impact our results given the different systematic effects (compared to earlier experiments) and improved sensitivity at higher energies; further AMS-02 observations might also admit improvements to our positron and antiproton constraints (and changes in the relative sensitivity to multi-step cascades) via improvements in the background and signal modeling.
10 An example application: the galactic center excess As an example application, we use our present results to explore the indirect detection constraints on DM explanations for the Galactic Center Gamma-Ray excess observed by Fermi, in the context of multi-step cascade annihilations [93]. Table 1 summarizes the results of [93]: the best-fit values of DM mass, cross section and step number for annihilations to electron, muon, tau and b-quark final states (note that other final states considered in the present work will be similar). In figure 6 we show these best-fit points for each channel, and the corresponding 1, 2 and 3σ confidence contours, together with the indirect detection constraints for that step number and channel. For comparison, we also show the indirect detection constraints for the case of direct annihilation; in this case the best fit occurs for annihilation to b-quarks [32]. (The best-fit points for direct annihilation to τ + τ − and µ + µ − final states occur at m χ ≈ 10 GeV and 5 GeV respectively, with corresponding cross-sections of 3 × 10 −27 cm 3 /s and 1.6 × 10 −26 cm 3 /s.) In all cases, the best-fit points lie very close to the bound from the Fermi dwarf observations. While the constraints can shift as a function of step number, the best-fit cross section and mass shifts in such a way as to maintain this borderline tension. This is not surprising, as Fermi observations of the Galactic center and dwarf galaxies probe exactly the same signal, i.e. photons within the energy range of Fermi. Any model that reproduces the GeV excess will yield a similar signal in dwarf galaxies.
The Planck and AMS-02 bounds, in contrast, do not share this degeneracy. However, for the e + e − and µ + µ − states where these constraints dominate, the region of interest for the Galactic Center excess is strongly excluded by both the robust Planck limits and the AMS-02 positron bounds. Figure 6. The yellow dot corresponds to the overall best fit point to the Galactic Center Excess for a given channel; electrons, muons, taus, and b-quarks computed in [93]. For best fit point we also display the 1,2 and 3 σ bands in gray. We have overlaid this with (95% confidence) constraints on σv corresponding to the best fit step shown in table 1. We also show for each channel the constraints for the case of direct annihilation to DM, and for the τ and b-channels we plot (orange point) the overall best fit for direct annihilation corresponding to m χ ≃ 46.6GeV and σv ≃ 1.60 × 10 −26 cm 3 sec −1 for direct annihilations to b-quarks and m χ ≃ 9.96GeV and σv ≃ 0.337 × 10 −26 cm 3 sec −1 for direct annihilations to τ + τ − [32].

JCAP06(2016)024
Final State n-step m χ (GeV) σv (cm 3   In order for indirect constraints on DM interpretations of the Galactic Center GeV excess to weaken significantly with increased step number, we see that the dominant bound should not arise from Fermi observations of dwarf galaxies. This could occur, for example, in the case where the final state is a 1 : 1 : 1 admixture of e + e − , µ + µ − and τ + τ − , with the τ + τ − component producing the gamma rays required to fit the GeV excess. Such a scenario has been proposed to explain hints of a low-energy counterpart to the excess, arising from inverse Compton scattering and bremsstrahlung of the electrons and positrons [94]. In the case of direct annihilation, limits from AMS-02 on the e + e − component severely constrain this scenario.
However, if we instead consider a 2-step cascade, which provides the best fit to the GeV excess for a τ + τ − final state, we should instead examine the Planck and AMS-02 bounds on 2-step annihilation to muon and electron final states with a DM mass of ∼ 24 GeV, with a cross section of ∼ 1.4 × 10 −26 cm 3 /s. None of the bounds we have derived on the two-step cascade can clearly exclude this scenario, as shown in figure 7.

Conclusion
We have shown that results from current DM indirect searches can be extended to constrain a broad space of dark sector models. We summarize our main points below: • Photon rich final states are generally most constrained by bounds from the Fermi dwarfs.
• Electron and muon final states are generally most constrained by AMS-02.
• The CMB bounds from Planck are robust and insensitive to number of dark sector steps. As a result the CMB bounds may become the most limiting in certain cases where the -21 -

JCAP06(2016)024
AMS-02 or dwarf bounds weaken as a result of large m χ or increasing number of dark sector steps.
• We find that for a fixed DM mass and final state, the presence of a hidden sector can change the overall cross section constraints by up to an order of magnitude in either direction (although the effect can be much smaller).
For hadronic SM final states (b-quarks, gluons, gauge bosons, Higgses), constraints from gamma-ray studies of dwarf galaxies generally remain the most limiting, and -within the kinematically allowed region -are generally fairly robust, although they can weaken at low masses and strengthen at high masses. More specifically for small but kinematically allowed masses the bound for final state gauge bosons, Higgses and b-quarks can weaken by about 0.1 decades. For the gluon final state, where very low DM masses are in principle possible, this bound can weaken by up to 1.1 decade; however a careful consideration of this regime would require taking into account the mass of the mediators, which may be comparable to Λ QCD . At high masses the bounds will strengthen by about 0.3-0.5 decades for the hadronic final states. The photon-rich tau and photon final states behave similarly, with the dwarf limits dominating the constraints except perhaps at very high masses (where it may be important to take constraints from VERITAS and H.E.S.S. into account). Adding extra cascade steps has little effect on the dwarf constraints on the photon final state at low masses (after the addition of the first cascade step, which weakens the limit by up to 0.8 decades), whereas for the tau final state it can weaken the bound by about 0.1 decades within the kinematically allowed regime.
For leptonic final states with few photons (electrons and muons), constraints from AMS-02 often appear to dominate the limits, but are quite sensitive to the number of cascade steps (as well as assumptions on the cosmic-ray propagation and local magnetic field; our limits are more conservative than others in the literature). At low masses (below a few hundred GeV), increasing the number of cascade steps can weaken the constraints by up to 2 orders of magnitude, at which point bounds from the CMB become more constraining. Above a few hundred GeV, however, adding more cascade steps tends to strengthen the constraints, so using results quoted for direct annihilation gives conservative bounds; the CMB limits are also generically stronger than the AMS-02 limits in this mass range.
If a quick estimate of constraints is needed, the CMB limits almost always appear to be within an order of magnitude of the strongest limit, for the cases we have tested, and vary by at most a factor of 1.3 to 1.5 over cascades with up to 6 steps.
The details of our code for n-step cascades which were used to produced our results are described in appendix D and are available at: http://web.mit.edu/lns/research/CascadeSpectra.html.

A Details of n-body cascades
In this appendix we will derive eq. (4.1) and provide some additional intuition for this case as well as pointing out that for a small number of steps, the cascade setup can provide an excellent approximation (albeit with some dependence on the channel). To set up this problem, firstly recall that the key physics encapsulated in eq. (2.3) is that when we add in a cascade step we need to boost the spectrum to the new rest frame. In the case of 2-body decays this is particularly simple, because we know exactly how much to boost by. Explicitly, if we have added in a step of the form φ i → φ i−1 φ i−1 , then in the φ i rest frame we know the φ i−1 particles must be emitted back to back, meaning we know their energy and hence their boost. If instead we introduce a step via φ i → φ i−1 φ i−1 φ i−1 , we no longer know the boost Figure 14. Dwarf limits for n-body vs m-step cascades for the γγ final state. We show the multibody case for n = 2-10 in gray, with lighter gray corresponding to larger n. In orange, green and purple we also show the 1, 3 and 5-step 2-body cascade for the same final state. As discussed in the text, for the multi-body case the spectrum sits in between the cascade spectra, and thus we expect the limits to do the same. The figure makes this clear and emphasises how the multi-body framework is captured within the cascade setup.
exactly, instead we can only associate a probability with any boost which we can determine from the energy distribution for a given φ i−1 . Accordingly what we need to calculate is the energy spectrum of a particular φ in the decay χχ → n × φ, and then combine this with a version of eq. (2.3) suitable for a general boost. Below we will firstly do this exactly for the case of a 3-body decay, show what this becomes after applying the large hierarchies approximation, and then we will show the general n-body result assuming hierarchical decays. As discussed, our starting point is the energy spectrum of a particular φ in the decay χχ → 3 × φ, which can be determined from the three body phase space. For this purpose we make use of the analytic formula for the n-body phase space outlined in [95,96]. In the case where our three final state scalars have mass m, we can write the 3-body phase space as: where λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2yz − 2zx and if we say the mass of the DM is m χ and the energy of one φ particle is E, then M 2 3 = 4m 2 χ and M 2 2 = 4m 2 χ + m 2 − 4m χ E. Using this, the energy spectrum of the scalars is simply: where the constant of proportionality can be determined by normalising the spectrum. Before proceeding, it is useful to introduce a set of dimensionless variables to work with as we did in the 2-body case. As there, we firstly define ǫ 1 = m/m χ , but note here that ǫ 1 ∈ [0, 2/3], rather than [0, 1] as in the 2-body case. To play a similar role to x, we also introduce ξ = E/m χ ∈ [ǫ 1 , 1 − 3ǫ 2 1 /4], where the limits here are fixed by eq. (A.1) and can also be seen from the kinematics. In terms of these variables, we can use eq. (A.2) and eq. (A.1) to arrive at:

JCAP06(2016)024
Figure 15. f eff for n = 0-6 step cascades to final state electrons and muons, with ǫ f = 0.3. Note that the difference in pattern between f eff for direct and single step electrons and higher step cascades can be understood by recalling that the direct electron FSR spectrum is sharply peaked. Each subsequent cascade step smooths out this spectrum thus changing the shape significantly.
where C is a constant that normalises the spectrum and can be determined numerically. Note that when ǫ 1 → 2/3, this distribution approaches a δ function, as expected when the particles are all produced at rest. We will return to the limit of small ǫ 1 shortly. Using this result, we can then revisit the derivation of the boost formula given in [13], a hierarchical version of which is given in eq. (2.3), and derive the analogue for an arbitrary boost. Doing so, if we label the spectrum of the decay of φ → 2 × (SM final state) as dN/dx 0 , we can write the spectrum of the same particle from the decay χχ → 3 × φ → 6 × (SM final state) as: where we have defined: There are two directions the above result can be generalised. For one, we could extend this to a longer cascade of 3-body decays, although the logic here is identical to the general 2-body case discussed in [13], so we will not repeat that here. Secondly we can look to extend this to 4-body decays and higher. The difficulty with this is that the n-body phase space quickly becomes analytically intractable. Nevertheless as observed in [20], in the large hierarchies regime (ǫ 1 ≪ 1) we regain analytic control as we will now outline. Returning to eq. (A.3), taking the ǫ 1 → 0 limit we find that: where now ξ ∈ [0, 1]. Following [20], this can then be generalised to the n-body case, where we find: where again ξ ∈ [0, 1]. Using this we can finally give the equivalent expression of eq. (2.3) for the n-body case: thereby demonstrating eq. (4.1). As a simple example of how this can be used, consider the decay φ → γγ which has the spectrum dN γ /dx 0 = 2δ(x 0 − 1). If we substitute this in, we find the spectrum for χχ → n × φ → 2n × γ is just Integrating this over x 1 ∈ [0, 1], we find N γ = 2n, as expected.
To follow on from this, consider the spectrum derived by repeated application of the boost formula in eq. (2.3) to the same φ → γγ spectrum, dN γ /dx 0 = 2δ(x 0 − 1). Doing so we obtain: Note that if we integrate this over x n ∈ [0, 1], we find N γ = 2 n+1 . Now by definition eq. (A.9) with n = 2 is identical eq. (A.10) with n = 1, as in this case they both represent a 2-body 1-step cascade. Note also though that if we take eq. (A.9) with n = 4 and eq. (A.10) with n = 2, then both situations have the same number of final state photons from different kinematic setups. In figure 12 we compare an n-step 2-body decay and a 1-step 2n-body decay for final state photons, for n = (1,2,3,4). We see that whilst they agree for n = 1 (by definition), and are quite similar to each other for n = 2, this similarity breaks down rapidly. This is not entirely surprising as a 6-body 1-step cascade has a different number of final state photons -30 -

JCAP06(2016)024
to a 2-body 3-step cascade, but one can also check that this latter spectrum also does not agree well with an 8-body 1-step result which would have the same number of photons. In figure 13 we show the same comparison for final state bb with ǫ b = 0.1. Here we see the agreement is better although still beginning to break down for larger n. We also tested this for several other final states, with the common theme that the spectrum of a 4-body 1-step decay is often well approximated by that of a 2-step 2-body decay. Lastly let us confirm the claim from the main text that the results for multi-body decays sit in between our multi-step cascade results. For this purpose consider again the photon spectrum obtained from φ decaying into bb with ǫ b = 0.1. In this case we plot the n-body spectrum for n = 2 to 10 in figure 3, which we presented in the main text. In the figure we also plot the case of a 1-step, 3-step and 5-step 2-body cascade. Clearly the n-body results sit in between these multi step cascade cases, which indicates that the constraints on an n-body 1-step cascade will be largely contained within the limits on a 2-body n-step cascade. To demonstrate this, in figure 14 we show that for the case of the γγ limits extracted from the dwarfs, the n = 2-10 multi-body decays sit exactly in between the 1, 3 and 5-step cascades as claimed.

B Details of the CMB results
In this appendix we present additional results from the CMB analysis. For many of the final states considered in the main text, the kinematic threshold on their production means that it is not sensible to go to lower masses than we presented. This is not the case, however, for electrons and muons, where we can take our results to much lower masses.
In figure 15 we present the value of f eff for cascades ending in final state electrons and muons with ǫ f = 0.3. Here we consider DM with mass as low O(keV) which is relevant to various CMB studies, which should be compared with our general results for f eff in figure 18 for DM annihilations into the eight final states considered in the main text. As expected f eff is largest for annihilations to final state electrons.
The corresponding bound σv for a given m χ for light DM annihilating to final state electrons and muons is displayed in figure 16, and more generally in figure 8. As this bound is fairly insensitive to the final state and number of steps it is interesting to examine the rescaled bound on σv /m χ which we display in figure 19 for n = 0-6 step cascades to various final states. We find that generically the re-scaled bound for all SM final states, ǫ f values, and cascade step falls within the narrow range σv /m χ = 10 −27.3 − 10 −26.6 cm 3 /s/GeV.

C Pass 7 versus Pass 8 for the dwarfs
As discussed in the main body, the limits displayed in figure 9 were derived using 6 years of Pass 8 data collected using the Fermi Gamma-Ray Space Telescope; more specifically using the publicly available results of [37] made from analysing this data. This work was an updated version of the analysis that appeared in [64], which set limits using the same 15 dwarf spheroidal galaxies, but only with 4 years of Pass 7 data. These results are also publicly available, 9 meaning we can cross check how much our results change when between datasets. We did this for each of the final states considered in figure 9 and found generically the shape of the limit curves were unchanged, but that the limits themselves improved by roughly half an order of magnitude when using the updated analysis. We show an example 9 Both can be obtained from http://www-glast.stanford.edu/pub data/.  figure 9 for the case of final state e + e − , for the case of the 4 years of Pass 7 data analysed in [64] (left) and for the 6 years of Pass 8 data considered in [37] (right). We see that the updated dataset essentially just strengthens the limits by roughly half an order of magnitude, without noticeably changing other basic features.
of this for the case of electrons in figure 17, and we see that the generic features of the limits are unchanged but the results strengthen as we move from the Pass 7 to the Pass 8 dataset.

D Description of cascade spectra files
All of the spectra used in this work are publicly available in .dat format at: http://web.mit.edu/lns/research/CascadeSpectra.html. The details of how these spectra were generated from the direct spectra mentioned in section 3 is discussed in section 2 and more comprehensively in [13]. The format of the spectra files has been modeled after those made available by [48], in the hope that anyone who has used the results of that paper should have no difficulty using ours. In addition to the files themselves we have also included two example files showing how to load the spectra in Mathematica and Python. There are four basic file types included, which we describe briefly in turn.
• AtProduction {gammas,positrons,antiprotons}.dat: these are the files provided by [48] and contain the 0-step or direct annihilation spectrum of {photons, positrons, antipro-tons} for various final states; • Again we emphasise that the AtProduction files were created by the authors of [48], we only include them in our results as it is convenient to store the 0-step spectra in separate files from the cascade results, yet having them in the same place is useful. As for the contents of the files, firstly the three AtProduction {gammas,positrons,antiprotons}.dat have the following format:   • Each file has 30 columns and 11099 rows, where the first row contains column labels and all others contain numerical values.
• The first column contains the DM mass in GeV, running from 5 GeV up to 100 TeV. Note using these direct spectra below 5 GeV is not advised as the extrapolation is often unreliable.
• The second column contains log 10 (x) values, where x = E/m χ . This ranges from -8.9 to 0 in steps of 0.05. Within these parameter ranges the interpolation is quite reliable, but outside these ranges linear interpolation is recommended. Note that several spectra, such as the γγ photons spectrum or the electron positron spectrum have no dependence on ǫ f or m φ . Nevertheless we still include an ǫ f column in those files for consistency, and note picking any value of this parameter will result in an identical spectrum.
• The second column contains log 10 (x) values, where x = E/m χ . This ranges from -8.9 to 0 in steps of 0.05.
• Finally the columns 3-8 contain the value of the spectrum in dN/d log 10 (x) = ln(10)xdN/dx of the spectrum at that value of m χ and ǫ f or m φ . The columns represent an n = 1 cascade (column 3) up to an n = 6 one (column 8).