On the Cosmic-Ray Spectra of Three-Body Lepton-Flavor-Violating Dark Matter Decays

We consider possible leptonic three-body decays of spin-1/2, charge-asymmetric dark matter. Assuming a general Dirac structure for the four-fermion contact interactions of interest, we study the cosmic-ray electron and positron spectra and show that good fits to the current data can be obtained for both charged-lepton-flavor-conserving and flavor-violating decay channels. We find that different choices for the Dirac structure of the underlying decay operator can be significantly compensated by different choices for the dark matter mass and lifetime. The decay modes we consider provide differing predictions for the cosmic-ray positron fraction at energies higher than those currently probed at the PAMELA experiment; these predictions might be tested at cosmic-ray detectors like AMS-02.


I. INTRODUCTION
Cosmic rays have been studied extensively at various earth-, balloon-and satellite-based experiments. Recently, the PAMELA satellite has observed an unexpected rise in the cosmicray positron fraction from approximately 7 to 100 GeV [1]. This feature is not explained by the expected background from the secondary production of cosmic-rays positrons. Moreover, observations of the total flux of electrons and positrons by Fermi-LAT [2] and H.E.S.S. [3] also show an excess over the predicted background, up to an energy of ∼ 1 TeV. The presence of nearby pulsars could provide an astrophysical explanation for these observations [4,5].
Nevertheless, more exotic scenarios remain possible. The annihilation of dark matter in the galactic halo to electrons and positrons provides one such possibility, though generic annihilation cross sections must be enhanced by a large boost factor in order to describe the data [6,7]. Alternatively, the excess could be explained by a TeV-scale decaying dark matter candidate. (For example, see Ref. [8]; for a recent review, see Ref. [9].) In this scenario, fits to the cosmic-ray data indicate that the dark matter must decay primarily to leptons with a lifetime of O(10 26 ) s.
While the thermal freeze-out of weakly-interacting, electroweak-scale dark matter can naturally lead to the desired relic density, this is not the only possible framework that can account for the present dark matter abundance. Recently proposed asymmetric dark matter models relate the baryon or lepton number densities to the dark matter number density, motivated by the fact that these quantities are not wildly dissimilar [10][11][12][13]. TeV-scale asymmetric dark matter models have been constructed, for example, in Refs. [11][12][13]. The asymmetry between dark matter particles and antiparticles can lead to differences in the primary cosmic-ray spectra of electrons and positrons, with potentially measurable consequences [14,15]. Evidence for such charge asymmetric dark matter decays would disfavor the pulsar explanation of the e ± excess [15]. In addition, charge asymmetric dark matter decays may allow one to discern whether dark matter decays are lepton-flavor-violating [16].
For example, the cosmic-ray spectra that one expects if dark matter decays symmetrically to e + µ − and e − µ + are indistinguishable from those obtained by assuming flavor-conserving decays to e + e − and µ + µ − with equal branching fraction; the same is not true if the dark matter decays asymmetrically to e + µ − alone, 100% of the time.
Refs. [15] and [16] study the cosmic-ray e ± spectra assuming a number of two-body charge-asymmetric dark matter decays, with the latter work focusing on lepton-flavor-violating modes. In this paper, we extend this body of work to charge-asymmetric three-body decays and, in particular, to modes that violate lepton flavor. We assume a spin-1/2 dark matter candidate that decays via four-fermion contact interactions to two charged leptons and a light, stable neutral particle. For the present purposes, the latter could either be a standard model neutrino or a lighter dark matter component. Four-fermion interactions have a long history in the development of the weak interactions, and one can easily imagine that dark matter decays could be the consequence of operators of this form, generated by higher-scale physics. Moreover, the possible presence of a neutrino in the primary decay may lead to interesting signals at neutrino telescopes [17]. Unlike the two-body decays already considered in the literature, the precise energy distribution of the decay products is affected by the Dirac matrix structure of these contact interactions, which is not known (unless a model is specified). By considering the most general possibilities, we show that different choices for the Dirac structure of the decay operators defined in Sec. II can be substantially compensated by different choices for the dark matter mass m ψ and lifetime τ ψ ; while the best fit values of these parameters change, the predicted spectra are not dramatically altered. On the other hand, we find that the flavor structure of the decay operator has a more significant effect.
Assuming various lepton-flavor-conserving and flavor-violating decay modes, we compute the resulting cosmic-ray spectra, performing χ 2 fits to the data to determine the optimal dark matter masses and lifetimes. Like Refs. [15,16], we obtain predictions for these spectra at e ± energies that are higher than those than can be probed accurately now. Future data from experiments like AMS-02 [18] may provide the opportunity to test these predictions, and evaluate them relative to other interpretations of the cosmic-ray positron excess.
This letter is organized as follows. In the next section, we discuss the assumed form of the dark matter operators. In Sec. III, we present the results of our numerical analysis and in Sec. IV, we discuss our results and directions for future work.

II. FOUR-FERMION OPERATORS
We consider a spin-1/2 dark matter candidate ψ that decays to ℓ + i ℓ − j ν where i and j are generation indices and ν represents a light, neutral particle. We assume that ν is either a standard model neutrino or a secondary dark matter component that is much lighter than ψ and contributes negligibly to the relic density. In the present analysis, the exact nature of the light neutral state will be irrelevant since its effect on our results will come solely from kinematics. We focus on the simplest scenario, in which there are no additional decay channels involving the charge conjugate of ν, and consider the possible four-fermion operators that contribute to the decays of interest. We work directly with the operators that may appear after the standard model electroweak gauge symmetry is spontaneously broken; for any operator found to have phenomenologically desirable properties, one may easily construct a gauge-invariant origin after the fact. Note that the production of a neutrino in the primary decay may have interesting phenomenological consequences (see, for example, Ref. [17]), which provides a separate motivation for our three-fermion final state. Once this choice is made, the dark matter spin must be 1/2 if the underlying theory is renormalizable 1 .
The problem of parametrizing an unknown decay amplitude of one spin-1/2 particle to three distinct spin-1/2 decay products was encountered in the study of muon decay, before the standard model was well established. The most general decay amplitude M can be parametrized by [20] where p ± and p 0 are the momenta of the decay products, labeled according to their electric charge, and the O i , i = 1 · · · 5 are elements of the set of linearly independent matrices The c i and c ′ i are complex coefficients. Terms involving the contraction of spinor indices that link different pairs of spinor wave functions can be recast in the form of Eq. (2.1) via Fierz transformations. Since the final state particles are much lighter than the dark matter candidate (which is at the TeV scale), we can safely neglect their masses.
Since the neutral final state particle is stable, the energy spectra of electrons and positrons that are observed at cosmic-ray observatories are determined by the energy spectra of the the charged leptons, ℓ + and ℓ − , that are produced in the primary decay; this follows from the differential decay distribution where |M| 2 is the spin-summed/averaged squared amplitude. We evaluate this quantity exactly from Eq. (2.1) using FeynCalc [21], and compute the ℓ ± energy distribution by integrating over the neutral lepton energy E 0 . We find that the result contains terms quadratic and cubic in E ± ; however, since the distribution must be normalized to unity, the result has the following simple parametrization: The requirement that this expression remains positive over the kinematically accessible range 0 ≤ E ± ≤ m ψ /2 restricts the parameters ξ + and ξ − to fall within the range The ξ ± are generally complicated functions of the operator coefficients c i and c ′ i ; we provide these in the appendix. In the present analysis, however, the exact relations are not particularly important; by leaving m ψ and τ ψ as fitting parameters, one obtains very similar predicted spectra, independent of the choice of the ξ ± . The fact that some solution exists for any desired Dirac structure of the underlying four-fermion operator makes it potentially easier to construct explicit models. Though we reserve the task of model-building to future work, it is worth noting, for example, that the operator We computed the electron and positron spectra using PYTHIA [22], taking into account the energy distributions of the primary leptons ℓ + and ℓ − .

III. COSMIC-RAY SPECTRA
To compute the relevant cosmic-ray fluxes, one must take into account that electrons and positrons produced in dark matter decays must propagate through the galaxy before reaching earth. While modeling this propagation is now standard in the literature on decaying dark matter scenarios, we briefly summarize our approach so that our discussion is self contained and our assumptions are manifest.

A. Cosmic-Ray Propagation
Let r be a position with respect to the center of the Milky Way Galaxy. We assume the spherically symmetric Navarro-Frenk-White dark matter halo density profile [23] ρ where ρ o ≃ 0.26 GeV/cm 3 and r c ≃ 20 kpc. The production rate of electrons/positrons per unit energy and per unit volume is then given by where m ψ and τ ψ are the dark matter mass and lifetime, respectively, and dN e ± /dE is the energy spectrum of electrons/positrons produced in the dark matter decay. Let f e ± (E, r) be the number density of electrons/positrons per unit energy. Then, f e ± (E, r) satisfies the transport equation [24] We assume the MED propagation model described in Ref. [25] for which where ǫ = E/(1 GeV). The diffusion zone is approximated as a cylinder with half-height L = 4 kpc and radius R = 20 kpc. We require f e ± (E, r) to vanish at the boundary of this zone. The solution at the heliospheric boundary is then given by [26] f e ± (E) = 1 The Green's function, G e ± (E, E ′ ), can be found in Ref. [26]. The interstellar flux of electrons/positrons created in dark matter decays is then given by where c is the speed of light. where, as before, ǫ = E/(1 GeV).
At the top of the Earth's atmosphere, these fluxes must be corrected to account for the effects of solar modulation [28]. The flux at the top of the atmosphere (TOA) is related to the interstellar (IS) flux by where E IS = E TOA + |e|φ F and |e|φ F = 550 MeV.
The total electron-positron flux is given by where k is a free parameter which determines the normalization of the background electron flux. In our numerical analysis, we find that the best fit values of k never deviate by more that two percent from 0.84 and that fixing k at this value has a negligible effect on the goodness of fits and our predicted spectra. Therefore, we set k = 0.84 henceforth to reproduce the cosmic-ray spectra at low energies. The positron fraction is given by (3.12) In the propagation model described above, the only remaining undetermined quantities are m ψ , τ ψ , dN e + /dE and dN e − /dE. The electron and positron energy spectra, dN e + /dE and dN e − /dE, are determined by m ψ and by a set of parameters which we describe in the following paragraph.
We consider dark matter decays of the form ψ → ℓ + i ℓ − j ν where ℓ ± i is a charged lepton of the i th generation. There are nine such decay channels, and we require where the B(ℓ + i ℓ − j ν) are branching fractions. For decays involving more than one channel, where (dN e ± /dE) ij is the electron/positron energy spectrum for ψ → ℓ + i ℓ − j ν. In Sec. II, we showed that the energy spectra of the charged leptons in the decay ψ → ℓ + i ℓ − j ν are characterized by the ordered pair (ξ + , ξ − ), where 0 ≤ ξ ± ≤ 96. We also showed that (dN e ± /dE) ij is entirely determined by m ψ and (ξ + , ξ − ). For decays involving more than one decay channel (e.g., ψ → e + µ − ν and ψ → µ + τ − ν), we assume a constant (ξ + , ξ − ). Then, since the branching fractions are subject to Eq. Leaving m ψ and τ ψ as free variables, we find that our results are relatively insensitive to the choice of (ξ + , ξ − ). This is demonstrated for the pure decay ψ → τ + τ − ν in Fig. 1 where we show the envelope of possible cosmic-ray spectra; that is, when we sample the (ξ + , ξ − ) parameter space, we find that all of the resulting curves fall between those plotted in Fig. 1. For the example shown, m ψ varies between 6.5 and 8.5 TeV while τ ψ varies between 0.5 × 10 26 s and 0.7 × 10 26 s; the χ 2 per degree of freedom (χ 2 /d.o.f.) remains between 0.5 and 0.6. We performed the same analysis on the other decay scenarios discussed below and found a similar behavior. As such, we take (ξ + , ξ − ) = (48, 48) for the remaining results that we present.
Note that, for fixed m ψ and τ ψ , the total electron-positron flux -which does not distinguish between the two electric charges -is the same for any two decays belonging to the same class. For this reason, we require only one plot of the total flux in Fig. 3. We find that the χ 2 is relatively flat as a function of the branching fraction within each class of decays: over the range of possible branching fractions, we find that the χ 2 /d.o.f. varies by no more than ±10% from 1.2, 1.1 and 0.6, for ψ → e ± µ ∓ ν, ψ → e ± τ ∓ ν, and ψ → µ ± τ ∓ ν, respectively.
Different choices for the branching fraction within a given class describe the existing data well, but provide different predicted spectra that interpolate between the curves shown.
Note that the distinctive dip in the µ + e − ν and τ + e − ν positron fractions around 1 TeV is due to the hard electron produced in the initial decay; this greatly enhances the electron to positron ratio in the high energy bins, leading to a suppression in the positron fraction for fixed total flux.

IV. DISCUSSION
The results presented in the previous section show that a variety of possible leptonflavor-violating decay modes for a spin-1/2, charge asymmetric dark matter candidate can   describe existing data well, as quantified by the χ 2 per degree of freedom for the best fits to the data. Significantly, the results for the predicted positron fraction differ substantially for energies above ∼ 100 GeV, the maximum for which the PAMELA experiment is sensitive.
In some case, more precise measurement of the total electron-postron flux around 1 TeV may also provide a means of distinguishing these scenarios. Future data from experiments like AMS-02 [18], which can probe these energy ranges of the predicted spectra, may determine whether the possibilities discussed in this letter present viable descriptions of the cosmic-ray spectrum.
In the meantime, the present work suggests a number of directions for further study: In the case where the stable, neutral particle in the final state is a standard model neutrino, one could study whether the decays of asymmetric dark matter that we have considered could be probed at neutrino observatories like IceCube [17] . One could also study additional astrophysical bounds on the scenarios described, for example, from the extragalactic gamma ray flux [15]. One can also attempt to find preferred forms of the underlying four-fermion operators (whose effects were parametrized in the present analysis by ξ ± ) by studying the simplest and best-motivated models that provide for their origin. Work in these directions is in progress and will be described in a longer publication.