New energy spectra in neutrino and photon detectors to reveal hidden dark matter signals

Neutral particles capable of travelling cosmic distances from a source to detectors on Earth are limited to photons and neutrinos. Examination of the Dark Matter annihilation/decay spectra for these particles reveals the presence of continuum spectra (e.g. due to fragmentation and W or Z decay) and peaks (due to direct annihilations/decays). However, when one explores extensions of the Standard Model (BSM), unexplored spectra emerge that differ significantly from those of the Standard Model (SM) for both neutrinos and photons. In this paper, we argue for the inclusion of important spectra that include peaks as well as previously largely unexplored entities such as boxes and combinations of box, peak and continuum decay spectra.


Introduction
The search for Dark Matter (DM) by indirect detection is the subject of many studies.A large number of experiments have investigated the cosmic antiproton, positron, photon and neutrino spectra.Notable experiments include but are not limited to, AMS-02 [1], Fermi-LAT [2], Icecube [3], ANTARES [4], H.E.S.S. [5], Pierre Auger [6], and VERITAS [7] which have been measuring charged cosmic rays, gamma rays, and neutrinos for decades.With the upcoming construction of new experiments such as KM3Net [8], the CTA observatory [9], and GRAND [10] the sensitivity to potential neutral particles arising from DM annihilation or decay will increase significantly.Historically the investigated particle spectra resulting from annihilating or decaying DM have focussed on processes involving two Standard Model (SM) final state particles, e.g.DM DM → b b, DM DM → τ + τ − , etc, which subsequently undergo decay, fragmentation and hadronization to produce antiproton, positron, photon, and neutrino spectra.These spectra are produced by well-understood mechanisms, leading to detailed analyses that can produce strong limits on the properties of DM.Some notable examples are the potential AMS-02 antiproton excess and the Fermi LAT gamma-ray excess, both of which have been explained by DM annihilation with a DM mass in the O(100) GeV region [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27].However, if one considers possible extensions beyond the standard model (BSM), other, still largely unexplored spectra are possible, which may differ considerably from the known standard spectra.In this paper, we describe a new type of neutrino and photon spectrum in a largely model-independent way.It contains a combination of a well-defined peak, a box and a combination of neutrino/photon spectra produced by the decay, fragmentation and hadronization of SM particles.Of course, a clearly defined peak has been investigated by Icecube [28], ANTARES [29] and Fermi-LAT [30] among others, as it is easily obtained by DM particles directly annihilating or decaying into the relevant final state resulting in a comparatively clean signal.However, a box shape and, of course, the combination of all three features leads to a significantly different spectrum.These types of spectra have been largely overlooked in both experimental and phenomenological research.To facilitate the search for these spectra, we provide a code to obtain a user-defined non-standard neutrino or photon spectrum by specifying the appropriate parameters.This paper is structured as follows.First, we detail the physics of the non-standard spectra and provide the relevant expressions of the kinematics.Next, some example BSM models are given that can produce non-standard neutrino or photon spectra.Then we detail the working of the sampling code and its verification, after which we provide some elementary parameter sets that capture the most important features of the spectra in order to facilitate experimental searches.We end with our conclusions.

Theoretical background
There are only two types of neutral particles that can travel cosmic distances from a source to detectors on or near Earth, namely photons and neutrinos.In the following subsections, the kinematics of the spectrum of a particle will be discussed as model-independently as possible.Both neutral particles can be described by the same kinematics, since they are both massless or have negligible mass.The main differences, of course, lie in the possible models that can produce such spectra.However, no such assumptions will be made in the following subsections.Moreover, both DM annihilation and DM decay can produce cosmic rays.The kinematics for DM decay is identical to DM annihilation, with the only difference being that the initial energy for (non-relativistic) DM decay is M DM , while for DM annihilation it is 2M DM .Thus, to obtain the kinematic expressions for the DM decay from those for the DM annihilation, one simply has to replace M DM by 1  2 M DM In the following subsections, we assume the standard scenario that two DM particles annihilate to neutrinos in order to simplify the discussion.

The box
The simplest non-standard spectrum arises from two DM particles annihilating into two BSM particles, X, which subsequently decay into neutrinos: where ν is any SM neutrino.More complicated decays of X are of course possible, which will be discussed in section 2.3.Depending on the model, of course, the XX-pair can also be a XX-pair, and any ν may well be an ν.This kind of DM annihilation leads to a 'box' shape of the neutrino spectrum.In the rest frame of the particle, X both neutrinos produced by the decay of X have a clearly defined energy of M X /2, where M X is the mass of X.However, the neutrinos need to be Lorentz boosted into the center-of-mass frame of the annihilating DM particles as they are now evaluated in the rest frame of the particle X, which of course differs from the center-of-mass frame.Assuming that DM has zero momentum at annihilation, using spherical coordinates in the rest frame of X, and choosing the motion of X in the z direction and thus boosted in the z direction, the energy of the neutrinos results in: Here M DM is the DM mass.The cos(θ) term is due to the spherical coordinates in the rest frame of X.To obtain a uniform distribution of points on a sphere, i.e. an isotropic decay, cos(θ) must be sampled uniformly between 0 and 1, as opposed to sampling θ uniformly between 0 and 2π1 .This uniform distribution in cos(θ) results in the neutrino having an equal probability of having any energy within the bounds of cos(θ) = −1 and +1, thereby resulting in a flat 'box' shape.However, a particle X will in general not decay isotropically if, for example, it is polarized, or has an asymmetric coupling to different helicities [31].The probabilty amplitude for the decay of X must be proportional to ⟨m ′ , S|R(θ)|m, S⟩, with |m, S⟩ and |m ′ , S⟩ the initial and final state respectively, where S is the total spin of X, m the angular momentum of X in its flight direction, m ′ the helicity diffence between the decay products of X, and R(θ) a rotation operator.The complete results for various spin and helicity configurations are provided in [31], which we have implemented here.The shapes of the boxes are given by: Here E X is the energy of X in the CM frame, C m ′ are model-dependent positive normalized coefficients, , that determine how X couples to the various polarizations.Note that for S = 1/2 the flat box is regained for either C 1/2 = C −1/2 , or when X has no preferred helicity.For a vector X the flat box is regained when X is unpolarized.

The peak
Another process is the annihilation of two DM particles into a neutrino and a BSM particle X: The neutrino, which comes directly from DM annihilation, forms a clearly defined peak in the spectrum with an energy of

Alternative decay modes for particle X
In most realistic BSM models, a possible particle X does not have a 100% branching ratio into νν, but can have multiple different decay modes, e.g.BR[X → νZ] = 0.5 and BR[X → W + e − ] = 0.5.The specific decay modes depend, of course, on the details and parameters of the chosen BSM model.Here, the assumption is explicitly made that X can only decay into SM particles.The SM particles 2 will undergo fragmentation, hadronization, and decay, thereby also producing neutrinos or radiating off photons.A labelling is made of the neutrinos that come from the various stages of DM annihilation into neutrinos: the neutrinos that form the peak are called primary neutrinos, those in the box are called secondary, and all neutrinos that come from the fragmentation and hadronization of an SM particle are called standard neutrinos: The other products of fragmentation, hadronization and decay, such as e ± , p, etc., are ignored here.Of course, the process DM DM → XX has no primary neutrinos.The energy range of the box depends, of course, on how the particle X is produced and what decay modes it has: for DMDM → XX .
The production mode of X changes the total energy of the particle X at creation and thus η, while the decay mode changes the energy of the secondary neutrino, giving the term M 2 SM .This results in the total width of the box being Notably, from this equation it can be seen that for DM DM → νX the energy of the peak is equal to the width of the box if M SM = 0.In general ).So for boxes with the same width, the peak is at the same energy, regardless of the energy of the box.Or, put another way, the position of the peak directly provides the width of the box up to the specific decay mode of X.
Of course, the mass of the particles DM and X can also be expressed by the initial and final energy of the box: For this expression, the decay mode of X must be fixed to know M 2 SM .The standard neutrinos produced by the hadronization and fragmentation of the SM particle can best be determined with programs such as Pythia [32] and are not analytically determinable.However, pre-computed spectra of DMDM → SM SM can be used to sample these neutrinos.For DMDM → SM SM the energy of the SM particle is simply E SM = M DM .For X → SM 1 SM 2 , i.e. two different SM particles, the energy of particle SM 1 and SM 2 is respectively given by: Thus, to obtain the correct neutrino spectrum of SM 1 , the spectrum of DMDM → SM 1 SM 1 is used, where M DM = E SM1 , and similarly for the neutrino spectrum of SM 2 .Since this neutrino spectrum of SM 1/2 is in the rest state of X, any neutrino sampled from this spectrum must be Lorentz boosted into the CM frame of annihilating DM particles, identical to the neutrinos coming directly from the decay of X.Additionally, when using the pre-computed spectra of DMDM → SM 1 SM 1 to determine the spectrum of a single particle SM 1 , its spectrum dN /dE is overestimated by a factor of 2 and thus needs to be compensated for.Furthermore, any correlations between the annihilation products of the pre-computed spectra are assumed to be negligible [33,34].
The average number of primary and secondary neutrinos per collision can be straightforwardly determined via counting: Table 1: The average number of primary (peak) and secondary (box) neutrinos per DM annihilation.
Note that when X decays into νν, the number of secondary neutrinos is doubled and the number of standard neutrinos is zero for that decay mode.The number of standard neutrinos is determined by the integral of their spectrum, dN/dE and therefore varies from case to case.

Model examples
The only requirement to obtain these spectra is a DM particle capable of either self-annihilation or decay into a mediating particle X that couples to neutrinos and/or photons.Of course, there are some limitations on the possible interactions and decay modes of the DM and X particles.For example, if X is a spin 1 2 particle, then the decay mode X → νν is forbidden by conservation of angular momentum, making a pure box shape impossible.Since the shape of the neutrino spectrum is determined by its kinematics, the spin of the DM or the X-particle is irrelevant when it comes to the position of the peak or the box; only the possible couplings are affected.Two examples of candidates for a particle X coupling to neutrinos are: Z ′ bosons, introduced in a U Lµ−Lτ (1) [35][36][37][38] gauge extensions of the SM and unstable heavy neutrinos such as those in an inverse seesaw mechanism [39][40][41].A heavy neutrino ν H and DM particle ϕ DM could, for example, have a non-standard muon neutron spectrum via the decay chain ϕ DM ϕ DM → ν µ ν H with ν H → ν µ Z, and the Z boson gives a continuous neutrino spectrum.The number of potential models can also be much larger; additional gauge groups, different neutrino mass mechanisms or Higgs mechanisms specific to neutrinos could all produce peak or box shapes.Of course, the number of different DM candidates is very large.These include particles that are added manually, such as complex scalar DM, or those that arise as a consequence of the theory itself, e.g.neutralinos in supersymmetry.Note that a single DM species annihilating with itself cannot have a neutrino spectrum consisting only of primary and secondary neutrinos, i.e.DMDM → νX with X → νν.The initial state has an even number of fermions, while the final state has an odd number of fermions, which is of course impossible. 3For photons, of course, there are no such constraints.It is possible that a process such as DMDM → ννX with X → νν, but here no peak is formed because the energy of the primary neutrinos and X is not uniquely determined by the two-body phase space.However, a spectrum with a peak and a box is easily obtained via the decay of fermionic DM.For example, a heavy neutrino decaying into a ν and a scalar, which subsequently decays into two ν's.In terms of models that might produce non-standard photon spectra, there are a plethora [42][43][44][45][46] of models proposed [47] to explain the 2015 750 GeV di-photon excess in ATLAS and CMS.More specifically [48] and [49] detail some concrete models regarding photon boxes.It should be noted that if X is electrically charged, as is very possible for a particle coupling to photons, the possible energy range of the spectrum is limited by constraints on the mass of X, e.g. by LEP searches.This is in contrast to a particle X that couples to neutrinos via the weak force, which can more easily evade experimental searches.Remarkably, X does not have to be electrically charged to couple to photons.The neutral pion π 0 , for example, has a decay mode into two photons.A non-standard photon spectrum could, for example, be made by a pion-like BSM particle Π 0 and a DM particle ϕ DM by ϕ DM ϕ DM → Π 0 Π 0 and Π 0 → γγ. 4 The sampling code

Sampling procedure
The sampling of the neutrinos is done by sampling a primary, secondary or standard neutrino according to the probability N prim /N , N sec /N and N stand /N , respectively, where N = N prim + N sec + N stand .Here N prim , N sec , and N stand refer to the number of primary, secondary and standard neutrinos, as defined by the Feynman diagrams in (6).The primary neutrino is always sampled at E peak since its energy is fixed and therefore has no distribution.The secondary neutrinos are either sampled uniformly, when the box is flat, or by rejection sampling if X has a polarization, with the box shapes given by Eq. ( 3).The limits of the box are given by Eq. ( 7).The standard neutrinos are sampled from precomputed spectra, and subsequently Lorentz boosted into the CM frame of the annihilating DM particles as described previously.We use the precomputed spectra from [50][51][52].In order to sample from these precomputed spectra, we numerically integrate them in order to obtain their cumulative distribution function.Subsequently, this cumulative distribution is used in order to perform inverse transform sampling.Notably, secondary photons can give rise to additional particles due to QED showering effects [32].We take these effects into account by using the aforementioned precomputed spectra when a secondary particle is a photon.

Verification of spectra
Figure 1: Comparison of the spectra between MadGraph/Pythia (orange) and the sampling code (blue).The upper left and lower two spectra were generated with simplified models, while the upper right spectrum is an example point from a full model [53].The lower two plots show the same model, but evaluated at different energies.
Figure 2: A comparison between the spectra generated by MadGraph/Pythia (orange) and the sampling code (blue) for two different spectra.The spectra of these plots are the same as the bottom two in 1, but plotted logarithmically.Note that the box now has a slope, which is purely a binning artifact.Since logarithmically equidistant bins are wider for higher values of E in absolute terms, a flat box in linearly equidistant bins will then have a slope when binned logarithmically due to the difference in entries per bin.
In order to verify the accuracy of this sampler, we perform cross-checks with multiple spectra computed with MadGraph5 v3.1.1 [54] and Pythia 8.309 across a range of masses and decay modes.We deem a simple visual inspection of the spectra to be sufficient to validate the fidelity of the sampler.In figure 1 four different example spectra are shown of which the process is indicated in the relevant plot.The top left and bottom two plots are generated through a simplified model in which X is a spin-0 or spin-1/2 particle respectively, while the spectrum from of top right plot is used from [53].The decay mode of X of the two spectra from the simplified models both have a branching ratio of 100% into a single decay mode, while the spectrum in the top right has a more complicated decay mode of BR , which shows that the sampler can indeed handle more complicated decay modes.In Figure 2, the two lower spectra of Figure 1 are plotted on a logarithmic scale, showing that the sampler at low energies performs accurately.There is some erratic behaviour at E ≈ 10 −4 GeV, but this is simply an artefact of a finite number of Monte Carlo samples.It should be noted that the normalization of the sampler and the spectra produced by MadGraph and Pythia are different, as the former only samples a spectrum while the latter fully computes an annihilation process.However, this difference can be easily remedied by counting the average number of neutrinos per event, which then directly provides the relative normalization factor.The spectra calculated by MadGraph and Pythia are the sum of 100,000 iterations, while the sampler's spectra are the result of 1,000,000 samples, which are then normalized.The normalization factor for an arbitrary spectrum can either be determined by counting the number of particles in a peak or box, which is directly related to the average number of particles per DM annihilation/decay as given by Tab. 1, or can be provided by the sampling code via the norm const() function.This function provides the normalization factor (N prim + N sec + N stand )/N samples .The sampled spectra are not normalized, such that they can easily be used for experimental searches, while the normalization is readily performed for phenomenological purposes.Remarkably, although the spectra shown are O(1-10)TeV, there is no inherent scale for these spectra, since both M DM and M X are a priori both unconstrained.Similarly, the detectability of any spectrum depends strongly on its annihilation cross-section ⟨σv⟩, which is of course model specific and which we do not comment on here.

Required input
The code can sample DMDM → ν/γX, DMDM → XX, DM → ν/γX, or DM → XX in which X can have any decay products and branching ratios into SM particles.The required input is as follows: • How the spectra is produced: DM decay (1) or DM annihilation (2) • The particle type: ν e , ν µ , ν τ , or γ.
• The DM mass M DM in GeV.
If no value is provided X is assumed to be unpolarized.
• The decay modes of X with the format being 'BR,daughter1,daughter2'.The total branching ratio must sum to 1.An optional parameter C m ′ − C −m ′ can be passed in order to fix the shape of the box according to Eq. ( 3).When no argument is given C m ′ − C −m ′ is assumed to be 0.
• The number of samples from the spectrum.
• The path to the csv file where the points are saved.When no input is given, no save is made.The keyword 'plot' or 'logplot' can be entered to plot the sampled data linearly or logarithmically.
Acceptable daughter particles of X decay are vl, e, mu, tau, h, z, w, ga, u, d, s, c, b, t, and g.The neutrino type is only important for the production of standard neutrinos; especially, the spectrum of tau neutrinos can differ significantly as compared to the spectra of electron and muon neutrinos when all other parameters are identical.The input can be provided manually or via an input file that is passed as an argument.
It is noteworthy that the masses of the DM and X particles are not constrained, while the interpolated spectra used for the daughter particles are tabulated up to 10,000 GeV, so the sampler cannot capture the spectra of the daughter particles with energy greater than 10,000 GeV.One disadvantage of these non-standard spectra is, of course, the enlargement of the parameter space.In the typical indirect search for DM, with a fixed channel, only the mass of the DM particle is important.In these non-standard spectra, however, the dimensionality of the parameter space is at least 2, namely the DM mass and the mass of the BSM particle X.

Conclusion
In this study, we have presented a novel approach for neutrino and photon spectra that is largely modelindependent.We have developed a user-specified Monte Carlo sampler to efficiently sample these spectra, which can be found in this github repository.While previous experimental searches have focused on neutrino and photon spectra with lines (peaks) and continuous spectra, more complex spectra have been largely overlooked.It is important to note that the spectra we have discussed are not constrained by mass, allowing for arbitrarily high or low energy ranges.Furthermore, the inclusion of a two-dimensional parameter space that includes the mass of dark matter and the mass of a BSM particle X, increases the complexity compared to the typical one-dimensional parameter space.This sampler may prove valuable not only for conventional high-energy astroparticle searches, but also for low-energy searches.