Non-Standard Neutrino Spectra From Annihilating Neutralino Dark Matter

Neutrino telescope experiments are rapidly becoming more competitive in indirect detection searches for dark matter. Neutrino signals arising from dark matter annihilations are typically assumed to originate from the hadronisation and decay of Standard Model particles. Here we showcase a supersymmetric model, the BLSSMIS, that can simultaneously obey current experimental limits while still providing a potentially observable non-standard neutrino spectrum from dark matter annihilation.


Introduction
One of the big unsolved mysteries in current-day physics is the exact nature of dark matter (DM), which as described by the Lambda-CDM model should make up 85% of the matter content of the universe [1].A standard hypothesis is that DM is comprised of Weakly Interacting Massive Particles (WIMPs), which for example can naturally be provided by extensions of the Standard Model (SM) in which supersymmetry (SUSY) is imposed and R-parity conservation is assumed.Collider, direct, and indirect-detection experiments have not yet found any conclusive evidence for the existence of WIMPs.While a number of different studies have found excesses that may indicate the presence of DM, for example the AMS-02 antiproton over proton ratio [2][3][4][5][6][7][8][9], and the Fermi-LAT gamma-ray data [10][11][12][13][14][15][16][17][18], these excesses remain small and subject to large uncertainties regarding production and propagation [19][20][21][22][23][24][25][26][27].An attractive messenger particle that circumvents these uncertainties is the neutrino, as these particles can propagate freely through the Galaxy, thus providing a potentially clean signal due to their low interaction rate.This low interaction rate is of course also challenging for detection purposes; it is only in recent years that cosmic-neutrino detection experiments have become potentially sensitive to neutrinos originating from DM annihilation [28][29][30].These spectra typically arise from SM final-state particles, e.g.b b, τ + τ − , or ν e,µ,τ ν e,µ,τ , in which ν e,µ,τ ν e,µ,τ provide the most stringent bounds due to the monochromatic neutrino line [31].There have been studies showing that, at least in the Minimally Supersymmetric Standel Model (MSSM), the neutrino spectra can deviate somewhat from only those resulting from SM final-state particles [32].The construction of KM3NeT is expected to increase this sensitivity, especially since neutrino signals arising from the Galactic Centre could be measured with unprecedented precision, thus providing an opportune way of investigating DM via neutrino physics.The SM, defined with only left-handed neutrinos, and its minimal supersymmetric extension, the MSSM, both lack a mechanism for generating neutrino masses.One example of a model that incorporate neutrino masses (and a DM sector) is the B-L-extended supersymmetric SM with inverse seesaw (BLSSMIS) [33][34][35].The inverse seesaw provides a mechanism for naturally light neutrinos by introducing right-handed neutrinos and an additional neutrino field with a small mass term.In the SM B-L (baryon number minus lepton number) is an accidental symmetry and is always conserved, as opposed to B and L symmetry separately, which are for example violated in sphaleron processes [36].This apparent accidental symmetry can be generated naturally from a U(1) B−L gauge group.Neutrino masses in the BLSSMIS are generated via an inverse-seesaw mechanism [37][38][39][40][41][42].This allows us to study the effects of a non-minimal neutrino mass model on the spectrum of cosmic neutrinos originating from DM annihilation.This paper is structured as follows.The relevant parameters, particle content, and other details of the BLSSMIS are discussed in section 2. The scanning methodology and the cuts that are imposed upon all data points is shown in section 3. Various DM observables are discussed in section 4, starting with the relic density, before moving to the current direct-detection limits and then discussing the LHC limits, and ending with a more in-depth study into the neutrino spectra of annihilating neutralinos.Lastly, in section 5, the concluding remarks are made.

The BLSSMIS model 2.1 Particle content
The total particle content of the BLSSMIS is largely the same as the MSSM, with the difference being a new Z ′ boson arising from the added U(1) B−L gauge group, and additional particles in the neutrino, sneutrino, higgs, and neutralino sectors.There are nine majorana neutrino mass eigenstates in this model due to the inverse seesaw mechanism, of which three are light SMlike neutrinos and six are heavy neutrinos.Naturally, this also extends the sneutrino sector, which contains eighteen real scalar fields.Two new higgs singlets are added to provide a mass to the Z ′ boson, so consequently three new higgs mass-eigenstates (two CP-even and one CPodd) appear additionally to those of the MSSM.The neutralinos are extended from four in the MSSM to seven neutralinos in the BLSSMIS; one new bino-like field and two higgsino-like fields.

Superpotential and Lagrangian
The increased particle content of the BLSSMIS leads to additional fields in the superpotential.The two new higgs singlet fields are denoted by η and η.Furthermore, to implement the inverse seesaw mechanism, six fermion singlet fields are added: three right-handed neutrinos, and three fields that are only charged under U(1) B−L , collectively called s 2 .To make the theory anomaly free, three additional singlet fields are required, denoted by s 1 , which cancel the anomaly introduced by adding the s 2 fields.The s 1 fields will be integrated out.Naturally, all fields are implemented as superfields since this model is supersymmetric.The resulting superpotential then follows as Here Q, L are the left-handed quark and lepton superfields respectively, Û, D, Ê, ν are the right-handed up-type quark, down-type quark, lepton, and neutrino superfields and ŝ2 is the superfield belonging to s 2 .Furthermore, Ĥu , Ĥd denote the up-type and down-type higgs superfields, and η, η are the aforementioned two new higgs singlet superfields.The Yukawa 3 × 3 matrices are indicated with Y i where the subscript i indicates the corresponding fields.We break SUSY by assuming Super Gravity (SUGRA) with unification at the Grand Unified Theory (GUT) scale, except for the gaugino breaking parameters.We define the GUT scale as the energy at which the coupling constants of the gauge groups all have a single unified value.
The breaking terms at the SUSY scale, whose value here is defined as the geometric mean of the two stop masses, are obtained via the evolution of the renormalisation group equations (RGEs).The resulting Lagrangian containing the soft-SUSY breaking terms then reads The product between two doublets is defined as a contraction with the antisymmetric tensor ε a b .The m 2 i are the sfermion and higgs mass terms, in which i again indicates the relevant field.The sfermion mass terms are matrices while the higgs mass terms are simply scalars.These scalar-breaking terms all have a common value m 0 at the GUT-scale.The T i are trilinear couplings, defined as Y i A 0 , in which Y i is the Yukawa coupling and A 0 a common GUT scale term.The gaugino mass parameters M 1 , M 2 , M 3 , M B L and M BB L are usually defined via a common M 1 2 at the GUT scale, but in this work each gaugino mass has its own GUT-scale parameter, i.e.M 1 , M 2 , M 3 , M B L for the B, W , G, and B B L respectively.The mixing term between B and B B L , M BB L , is set to zero at the GUT-scale.The splitting of the gaugino mass term grants greater freedom in the neutralino sector and is allowed in SUGRA models as the unification of the different gaugino mass parameters is not required, merely often implemented to reduce the number of parameters.In addition to the standard tan(β) parameter, which defines the value for the ratio of the two higgs vacuum expectation values (vevs) v u /v d , there is an additional tan(β ′ ) that represents the ratio of the two new vevs, v η /v η .Furthermore, we fix |µ| 2 , |µ η | 2 , Re(B µ ), and Re(B η ) using the tadpole equations that follow from the minimisation of the higgs scalar potentials.Additionally, the charges under U(1) B−L are fixed for all fields except the new higgs fields η and η, and the new neutrino fields ŝ1 and ŝ2 .We choose the charges of ŝ1 and ŝ2 to be 1/2 and -1/2 respectively under U(1) B−L .This choice fixes the U(1) B−L charge of η and η to be -1 and 1 respectively.A more detailed discussion on the charges of the superfields and the corresponding superpotentials can be found in appendix A. The foregoing implies that the coupling strength of the Z ′ boson to all fields is fixed by the theory and cannot be scaled manually.The only impact that the model parameters have on interactions involving the Z ′ particle is the mass of the Z ′ , and the mixing and mass of the mass eigenstates of the particles coupling to the Z ′ .While at the GUT scale the couplings are unified, mixing between the U(1) Y and U(1) B−L occurs at the SUSY scale.This results in the gauge couplings in the covariant derivative appearing in the dynamic terms of the Lagrangian becoming a 2×2 matrix.These mixing effects are a priori expected to have an impact on the electroweak sector, but this can be removed by redefining the two gauge fields.The fields before mixing, B ′ and B ′ B−L , can be redefined into B and B B−L such that the electroweak sector is left untouched.This is done by performing a simple rotation, i.e. where Here g 1 and g B−L are the U(1) Y and U(1) B−L coupling constants respectively.At the GUT scale, all coupling constants are unified, while the off-diagonal elements are zero.This completely fixes the new coupling constants and thus they cannot be tuned freely.
The mixing of the U(1) gauge groups impacts the M B L L term, but also the neutralino sector as it introduces an explicit mixing term between the B B L and the higgsino fields.In the basis of B W h d h u B B L h η h η , the mass matrix of the neutralinos is given by Here the upper-left 4×4 matrix can be recognised as the MSSM neutralino mixing matrix.When refering to the gaugino content of the mass eigenstates, in what follows, we will combine the gauge eigenstates h d and h u into the total higgsino component h, defined through The h η and h η will similarly be combined into h B L .The neutralino mixing matrix is diagonalized as N * M χ 0 N −1 = M D χ 0 .The neutralino mass eigenstates are given by χ 0 1...7 , which are mass ordered.The chargino sector is the same as in the MSSM, and thus shall not be discussed here.The mass matrix for the neutrinos can be read from the superpotential (1) and in a basis of ν L ν R s 2 it reads Here Since the µ S parameter breaks the B − L symmetry it can only receive a small non-zero value via high-scale radiative effects [38].Assuming µ S to be small automatically leads to three light neutrino mass eigenstates (ν l ) and six heavy (ν h 1,2 ) mass eigenstates.Furthermore, taking M ν , M X , and µ S to be m ν I, m X I and µ S I respectively, where I is the 3 × 3 identity matrix, we find the tree-level neutrino mass eigenstates as The mass eigenstates of the heavy neutrino have no suppression of the left-handed neutrino field, and thus couple without suppression to the higgs, W ± , and Z boson.Additionally, a coupling to the Z ′ via the B-L charge of the left and right-handed neutrinos is present.Both aforementioned couplings result in the ν h being an unstable particle.We shall refer to both ν h 1 and ν h 2 as ν h for the remainder of the text, as they are phenomenologically equivalent for our purposes due to their similar mass, decay width and decay modes.Thus of the ab initio three DM candidates of the B-L-SSM-IS, i.e. the heavy neutrino, sneutrino, and neutralino, the heavy neutrino can be discarded by this assumption as a DM candidate.However, while excluded as a DM candidate, the heavy neutrinos play a significant role in the neutrino spectrum of annihilating neutralinos for which B B L or h B L is the largest component.Furthermore, in the following we shall only look at neutralino dark matter, as sneutrino dark matter in the B-L-SSM-IS has been sufficiently investigated in previous works where it is shown that it is indeed a viable DM candidate [33,34].

Scanning methodology
The objective of our scanning procedure was to find solutions with novel phenomenology, and not to completely map out high-likelihood regions in the parameter space.In order to do this we implemented the following search strategy.Three different initial parameter ranges were used, which are given in table 1.We scanned these parameter spaces using a particle filter [17,[43][44][45][46][47] in which we applied hard cuts on the exclusion limits, which will be explained below, in order to increase the efficiency of the sampling.The scanning procedure was performed in three steps.As a first step, the entire parameter space was sampled randomly and uniformly, with the exception of Y ν , Y x , and µ S which are sampled logarithmically.The dataset that this yielded was used as the seed for the second step of the sampling procedure, in which a Gaussian particle filter was used to find regions of the parameter space that provide a DM relic density of 0.12 or less.Having found these viable regions, we manually identified the interesting regions, namely those with novel neutrino phenomenology, contained in this dataset and fed those as new input into a Gaussian particle filter in order to zoom in on those regions of interest.In general we decreased the gaussian width by 0.5 per iteration, but this decrease is occasionally manually adjusted in order to keep a good resolution.We stopped the scan after we found solutions with new phenomenology.Note that µ S is chosen to be small, in accordance with its nature as discussed in the previous section.Furthermore, we choose the positive branch of both sign(µ) and sign(µ η ).The three ranges denoted in table 1 are used to investigate the high-mass parameter space, and no significant difference was found between the three ranges, aside from the presence of more high-mass solutions when the parameter ranges are increased.In total O(10 7 ) points have been considered.Sarah 4.14.3 [48][49][50][51][52][53] is used as input for SPheno-4.0.4 [54,55], which is employed as Table 1: The initial parameter ranges for the three different scanning ranges.The parameters are defined and sampled at the GUT scale.the spectrum generator.A Universal Feynrules Output [56] file was manually written for MadGraph.1 Micromegas 5.2.1 [58][59][60][61], which also uses the Sarah-generated files as an input, is used to compute the DM relic density Ωh2 , spin-dependent and spin-independent DM cross sections for the proton and neutron, σ SD respectively, and the velocityweighted DM-annihilation cross section 〈σv〉, in addition to the active (co)-annihilation channels.The direct-detection limits are implemented using DDCalc 2.2.0 [62] using the most recent limits of Xenon [63,64], PICO [65], LUX [66,67], and PandaX [68,69].The spindependent and spin-independent cross sections are scaled with (Ωh 2 )/0.12 ≡ ξ in order to account for any DM underabundance.Similarly, the velocity weighted cross section is scaled with ξ 2 .The scaled values are used to interpret the direct-detection limits on Ωh 2 , and the indirect-detection limits on 〈σv〉 from the Fermi-LAT gamma-ray limits of the Milky-way dwarf spheroidal galaxies [70] and the limits from IceCube and ANTARES on the neutrino limits from dark matter annihilation [29,30], as these limits typically are produced with the assumption that only one single DM species exists.Furthermore, to obey limits from LHC searches, the σ(pp did not do so already. 2 The LHC exclusion limits using the aforementioned production cross sections are determined with Smodels 2.0.0 [72][73][74][75][76]. The cuts in table 2 have been imposed on all generated model points as presented in section 4. The mass of the lightest chargino should be larger than 103.5 GeV, as charginos below this mass are excluded by LEP.The higgs sector has been expanded with two new higgs singlets, thus its mass computation obtains additional terms, introducing additional uncertainties.To be conservative, the SM-like higgs boson is required to have a mass between 120 and 130 GeV.The maximum variation of 5 GeV on the higgs mass is expected to have little impact on the DM candidates that we will focus on in what follows, the B B L and h B L -like neutralinos.Furthermore, both the direct detection and LHC limits both need to be satisfied.Note that especially the cuts on the squark, slepton en gluino masses are more stringent than strictly Table 2: Constraints that are imposed on all model points.The DM direct detection and the LHC production cross sections (σ(pp → χ 0 2/3 χ ± 1 ) and σ(pp → χ + 1 χ − 1 )) are implemented using DDCalc and Smodels respectively.The constraints on particle masses are implemented as a hard cut.Note that constraints on the sfermions and gluino masses are stricter than those provided by LHC searches, but this is done as so to guarantee that the model is not excluded by searches for coloured sparticles.Similarly, the Z ′ -boson mass is cut more stringently than required.
needed, but are chosen such that they easily evade the LHC limits and are factored out of the phenomenology.We deem these cuts reasonable for our purposes due to the limited impact of coloured sparticles on the DM sector.
We stress here that in the subsequent plots the number density of points is in no way indicative of any statistical significance, but simply a result of sampling; it is always possible to increase the sampling in the surrounding parameter space near low-density regions.

Relic density
The relic density in the case of neutralino lightest stable particle (LSP) can be seen in figure 1.
The vast majority of W -like and h-like LSPs are pure states.We consider a neutralino to be a pure state when the total fraction of any component exceeds 0.99.The pure states lie on a line with an M 2  the g × coupling constant.The B B L -h mixing term can clearly be seen in the neutralino mass matrix of Eq. (5).A funnel is also present for B, B B L and h B L with the second lightest higgs h 2 .However, the funnel cannot be seen in figure 1 since the mass of the second higgs is variable.Similarly a funnel for the third CP-even higgs h 3 and first CP-odd higgs A 0 1 exists, but these funnels are again not visible due to their variable masses.A funnel for h 4 and A 0 2 could not be found due to the high mass of both particles, since no solutions were found in which the mass of the first neutralino is approximately half that of h 4 or A 0 2 .Moreover, many of the B-like LSPs that result in a relic density that does not exceed the observed value are those where the first chargino and second neutralino have similar mass to the first neutralino.This enables a χ ± 1 and χ 0 2 (co)-annihilation mechanism in the early universe, thereby sufficiently lowering the relic density [83].
Figure 2 shows the B B L (left) and h B L (right) components for all LSPs that have passed the constraints of table 2, including a constraint on the relic density, Ωh 2 < 0.15.Interestingly, of all B B L -like LSPs, none are pure B B L states.This is expected from the neutralino mass matrix Eq. ( 5), which contains explicit mixing terms between the B, h, and h B L fields.Similarly many h B L -like LSPs have a significant B B L mixture.This high degree of mixture is caused by the relation between the mass of the Z ′ boson and the mass of the LSP: the off-diagonal terms mixing the B B L and h B L fields in the neutralino mass matrix of Eq. ( 5) contain the terms −g B−L v η and g B−L v η , thus when the mass of a h B L neutralino gets close to the that of the Z ′ boson, these terms become relevant and mixing between B B L and h B L will occur.Of the found h B L -like LSPs the majority of the solutions with a relic density of 0.15 or lower have a funnel

Direct detection limits
Figure 3 shows the spin-dependent and spin-independent cross sections, relevant for DM direct detection, of all non-excluded points with Ωh 2 < 0.15. 4 The tentative LZ limits [84] are shown as a dotted line. 5It can be seen that the SM-like Higgs funnel will be excluded, and the higgsino and BBL region is probed.Additionally, the projected PICO limits for the spin-dependent cross section and the Xeno n-nT and Darwin limits for the spin-independent cross sections are shown.It can be seen that the projected limits on the spin-dependent cross section will not probe any of the found solutions.On the contrary, the projected spin-independent limits are expected to probe all of the found pure higgsino solutions.Similarly, the B B L -like LSP solutions will be within reach of Xenon-nT and DARWIN, but there remain solutions outside the reach of these experiments.It is therefore expected that, at least for the h B L and some B B L -like LSPs, planned direct-detection experiments will not provide conclusive evidence on their potential existence.For DM-nucleon scattering the value of the spin-dependent cross section is driven by a tchannel exchange of a Z or Z ′ boson, or a t-channel exchange of a squark.Similarly, the spin-independent cross section is driven by a mediating higgs, which can be any of the four CP-even higgs mass eigenstates, or again a mediating squark.Note, however, that the squark contribution to either the spin-dependent or spin-independent cross sections is often negligible due to their high masses of O(2-10 TeV) for the found solutions, as the squark masses are cut at 2 TeV.The squark contributions become relevant for small χ 0 1 χ 0 1 Z or χ 0 1 χ 0 1 h couplings.However, for these scenarios, both the spin-dependent and spin-independent cross sections are much lower than the current or projected limits, so these models are not relevant for nearfuture DM phenomenology, and we will not discuss them further.The higgsinos feature the largest spin-dependent cross section due to their coupling strength with the Z boson.However, this coupling is suppressed for pure higgsino states since it is proportional to N * 13 N 13 − N * 14 N 14 , which for pure states becomes equal to zero.Notably, this suppression is not present for the χ 0 1 χ 0 1 Z ′ coupling; while this coupling contains a similar term of the form N * 16 N 16 − N * 17 N 17 , the h η and h η components are in most cases not the same, even for pure h B L states, due to the mixing terms between the B B L and h B L fields in the neutralino mass matrix.The bino and wino-like LSPs have a smaller coupling strength to the Z boson compared to the higgsino LSPs, as these solutions by definition have a less higgsino sizable component compared to the higgsino neutralinos.The spin-dependent cross section for both the h B L and B B L -like LSPs is generated via a mediating Z and Z ′ .The Z ′ contribution to the spin-dependent cross section for both the B B L and h B L has two dependencies: the amount of h B L present in the neutralino and the mass of the Z ′ .Given the heavy Z ′ -boson mass, the spin-dependent cross section for the B B L and h B L -like LSPs is typically low.Similarly to the spin-dependent cross section, the higgsino solutions also feature the highest values for the spin-independent cross section of all LSP types.Naturally, this results from their relatively large higgsino-bino-higgs and higgsino-wino-higgs coupling.However, in contrast to the spin-dependent cross section, both the B B L and h B L LSPs generally have higher cross sections than the wino LSPs.The cross section for models with a B B L -like neutralinos is driven mainly by the SM-like higgs.The h B L neutralinos couple to the new η and η-like higgses as long as a sizeable amount of B B L is present, which is typically the case as seen from figure 2.
However, the size of the spin-independent cross section is suppressed by the masses of the new higgses, which can be heavy as they are driven by the Z ′ mass, which in turn needs to be heavier than 2.5 TeV.The wino and bino solutions mostly have a low spin-independent cross section due to the correction factor on the direct detection cross sections of any relic density underabundance ξ = Ωh 2 /0.12.

LHC production cross sections
While the LHC phenomenology of the BLSSMIS is not focus of study in this paper, we here briefly comment on the typical production cross sections that our spectra predict.Figure 4 shows the size of the pp → χ 0 2 χ ± 1 production cross section against the mass of the second neutralino, with the colour coding indicating the dominant contribution of the second neutralino.Unsurprisingly, two clear lines can be seen for a wino-and higgsino-like second neutralino, since the pp → χ 0 2 χ ± 1 cross section is the largest for wino-or higgsino-like second neutralino and first chargino.However, when both the first and second neutralino have no sizable wino or higgsino component the cross section is not correlated to the mass of the second neutralino, which indeed can be observed in figure 4. In the MSSM, the first or second neutralino is guaranteed to contain a sizeable wino or higgsino component or a mixture of both, which increases the chargino-neutralino production cross section.Moreover, in the MSSM the second neutralino cannot be a pure bino state, since the first neutralino then needs to be composed entirely of wino and higgsino components, which are excluded by either the relic density or dark matter direct detection experiments.This is no longer true in the BLSSMIS, as the B, B B L , and h B L components are able to compose the first four neutralino mass eigenstates, such that no large wino or higgsino component is Figure 4: The σ(pp → χ 0 2 χ ± 1 ) production cross section in pb against the mass of the second neutralino in GeV.The points are labelled according to the dominant contribution of the second neutralino.The shown points pass all constraints of table 2 and satisfy Ωh 2 < 0.15.
present for these neutralinos.This results in the possibility of much lower production cross sections as compared to the MSSM.While the h B L -like first LSPs that pass the relic density cut are typically heavy, B B L -like LSPs can be low-mass.From figure 2 one can infer that low-mass B B L -like LSPs have sizable B and h B L components.These scenarios typically translate to the lightest two neutralinos being composed of B, B B L , and h B L , resulting in a potentially very low pp → χ 0 2 χ ± 1 cross section.Thus the B-like second neutralinos in figure 4 correspond mostly to B B L -like LSPs.We have additionally computed both σ(pp → χ 0 3 χ ± 1 ) and σ(pp → χ + 1 χ − 1 ).Here σ(pp is taken into account due to the aforementioned possibility of the first two neutralinos to not have a sizeable wino or higgsino component, and thus for the third neutralino to be the first wino/higgsino neutralino.Wino/higgsino neutralinos typically feature higher production cross sections than bino-like neutralinos.We disregard the neutralino-chargino production cross sections that include the fourth, fifth, sixth and seventh neutralino, as these most likely will not exclude a given model point if it has not been done so by the previously mentioned cross sections due to a lower cross section, an increased complexity of the neutralino decay chain, or both. 6Since the focus of this paper lies on studying the phenomenology of our spectra at neutrino experiments, we refrain from commenting further on a possible LHC optimised search for these models.

Indirect-detection limits
From the previous subsections, we can conclude that most h B L -like DM solutions are not expected to be probed based on direct detection or collider experiments, and only future direct detection experiments may be sensitive to B B L DM.This leaves indirect detection as a very interesting, and maybe only possible, short-term viable detection method.The scaled velocity-weighted cross section can be seen in figure 5.The Fermi-LAT limits from dwarfspheroidal galaxies [70], HESS limits from the Galactic centre [78], and the IceCube and ANTARES [29,30] limits on the resulting neutrino spectra are implemented, and for specific annihilation channels indicated by the solid blue lines in figure 5.In order to be conservative, we do not include the limits on the velocity-weighted annihilation cross section arising from antiproton limits due to large inherent hadronization and propagation uncertainties [19][20][21][22][23][24][25][26][27].Furthermore, the projected KM3NeT limits for DMDM → νν is shown as a dotted line [85].It can be seen that KM3NeT will be able to reach both the χ0 1 χ0 1 → ν l ν h and νh νh regions.The size of 〈σv〉ξ 2 for B B L and h B L -like DM is largely determined by the presence of either a Z ′ or higgs funnel.If the corresponding DM masses are close to a Z ′ or higgs resonance the velocity-weighted cross section increases rapidly.The B, W , and h neutralinos annihilate predominantly into SM particles; W + W − , leptons, quarks and the SM-like higgs.The branching ratios of neutralinos annihilating into the different final states depends on the exact composition of the neutralino, but only the branching ratios of annihilating neutralinos into neutrinos, either light or heavy, is of interest here.The annihilation spectra of the bino, wino and higgsino neutralinos will not be investigated further as these neutralinos annihilate mostly into SM particles, and the interest here lies with neutrino phenomenology that is not present in the MSSM.Both the B B L and h B L -like DM can annihilate mostly (Br[ χ 0 1 χ 0 1 → νν] > 0.5) into neutrinos if either χ 0 1 χ 0 1 → ν h ν h or χ 0 1 χ 0 1 → ν h ν l is kinematically allowed.Both these processes are predominantly mediated via a Z ′ boson, 7 and thus B B L and h B L -like neutralinos 6 In the case that the third neutralino is the first higgsino-like one, the fourth neutralino will very likely be also higgsino-like.Resultingly, the production of pp → χ 0 4 χ ± 1 will be of a similar size to that of pp → χ 0 3 χ ± 1 .However, this will only increase the total production cross section by a factor of two, and is therefore only relevant for those spectra that feature production cross sections that are on the verge of detection.We deem these scenarios sufficiently unlikely that they can safely be neglected. 7Other mediating particles exist, but they are subleading in their contributions.can annihilate mostly into neutrinos.Again, the specific branching ratios depend on the exact composition of the neutralino.Notably, for all neutralinos that annihilate mostly into ν h ν l the process χ 0 1 χ 0 1 → ν h ν h is kinematically forbidden.Here ν l and ν h are the neutrino mass eigenstates as given in Eq. (7).In total two relevant DM annihilation channels exist when regarding the neutrino spectra: χ 0 1 χ 0 1 → ν h ν h , and χ 0 1 χ 0 1 → ν h ν l .In general the process χ 0 1 χ 0 1 → ν h ν l process where the heavy neutrino ν h subsequently decays is The momenta of the neutrinos resulting from neutralino annihilation can be computed using simple 2→2 kinematics.Moreover, the momenta of the neutrinos are completely fixed when assuming the velocity of the neutralinos to be negligible, which can safely be done seeing as the present-day DM particles are non-relativistic.For the process χ 0 1 χ 0 1 → ν h ν l the light neutrino ν l then has an energy of where M χ 0 1 is the mass of the neutralino and M ν h is the mass of the heavy neutrino.As the light neutrino has a single unique energy, rather than a distribution, a clearly defined peak in the neutrino spectrum must be present.The heavy neutrino, via its decay products, of course also impacts the neutrino spectrum.
In all cases studied the heavy neutrino has three main decay modes Br[ν h → ν l Z] ≃ Br[ν h → ν l h] = O(0.25)and Br[ν h → W ± l ∓ ] = O(0.5).The Z, h, W ± and l ± contribute to the neutrino spectrum via hadronisation and decays, and have a standard contribution to the neutrino spectrum.The contribution of SM particles to the total neutrino spectrum is best determined via Monte Carlo sampling.However, the contribution of the light neutrino can be computed analytically.In its rest frame, the heavy neutrino decays isotropically, thus the four-momenta of the two daughter particles are easily computed in general spherical coordinates in this frame.However, the heavy neutrino needs to be boosted from its rest frame to the rest frame of the annihilating neutralinos to determine the contribution of the light daughter neutrino to the total neutrino spectrum.When boosting in the z direction8 the energy of the daughter neutrino in the rest frame of the annihilating neutralinos is Here m h,Z denotes the mass of the Z or h boson featuring in the decay ν h → ν l Z/h.Note that the cos(θ ) term is from the momentum of the light neutrino arising from the decay of the heavy neutrino in the z direction in the rest frame of the heavy neutrino.Thus the energy of the light neutrino depends on the angle θ , which implies that the contribution of this neutrino to the total neutrino spectrum has an energy range that depends on θ .It should be noted that, since the general four momenta of the neutrinos are written in terms of spherical coordinates, the angle θ needs to be sampled according to cos −1 (1 − 2u) with u ∈ [0, 1] to get a uniform distribution of points on a sphere, 9 as is required by the isotropy of the decay of the heavy neutrino.This causes the energy spectrum of the light neutrinos coming from heavy neutrino decay to lie on a flat line, i.e. it will form a plateau.The edges of the plateau correspond to cos(θ ) = 0 for the high end and cos(θ ) = −1 for the low end as can be seen from Eq. ( 9).The number of neutrinos in the peak and plateau regions is given by Here N ν,peak and N ν,plateau are the number of neutrinos in the peak and plateau regions respectively.Furthermore, A and B are used to indicate any particle flavour (including heavy/light neutrinos) and N f the number of light neutrino flavours.Note that if A = ν h the number of neutrinos in the plateau is twice as large compared to the case where A ̸ = ν h .In figure 6 the neutrino spectra of two example spectra (whose LSP and heavy neutrino masses, plateau and peak energies are shown in table 3) are shown that have as the dominant annihilation channel χ 0 1 χ 0 1 → ν h ν h on the left and χ 0 1 χ 0 1 → ν h ν l on the right.The predicted peak and plateau regions are indicated by the shaded regions.Both features can clearly be seen in the computed spectra.The peak region can be seen at E ≈ 950 GeV for the ν h ν h channel and at E ≈ 300 GeV for ν h ν l .The plateau is visible between 300 ≲ E ≲ 850 GeV for ν h ν h and 500 ≲ E ≲ 800 GeV for ν l ν l .Note that the peak for the ν h ν h dominant channel in figure 6 is due to the presence of a sub-dominant annihilation channel χ 0 1 χ 0 1 → ν h ν l which has a branching ratio of approximately 1%.The computed values of the peak and plateau regions are shown in table 3. When compared to the spectra of figure 6 these correspond precisely to the MadGraph data,10 up to binning effects.Table 3: The numerical values of the energy of the peak (E peak ), the low-energy (E plateau low ) and high-energy (E plateau high ) boundary of the plateau of the neutrino spectrum for the two show-case files in GeV.Additionally, the height of the peak of the right (left) plot in figure 6 can be seen to be ∼0.3 (∼0.0026), which is as predicted by Eq. (10) since BR[ χ 0 1 χ 0 1 → ν h ν l ] is ∼0.96 (∼0.01) and BR[ν h → ν l B] ≈ 0.5.Similarly, the total number of neutrinos in the plateau of the right (left) plot is ∼0.16 (∼0.36), which is again as expected based on Eq. ( 11).From these two example spectra, it can be verified that the shape and size of the peak and plateau regions can accurately be predicted given the model parameters.Moreover, the shape of the peak and plateau regions is uniquely specified.Thus any possible future measurement will directly provide insight into the neutralino mass, the heavy neutrino mass, the branching fractions of a heavy neutrino, and the branching fraction of two annihilating neutralinos into ν h ν h and ν h ν l .We however, refrain from giving exclusion lines on this spectrum.While the peak of the spectrum can be directly tied to the monochromatic lines of ν µ ν µ exclusion limits, namely the exclusion limit for χ0 1 χ0

Main channel m χ
1 → ν l ν h is half that of DM DM → ν e,µ,τ ν e,µ,τ at the worst, the impact of the plateau and remaining features of the spectrum in conjunction with each other requires experiment-specific analysis, especially the detector response.Especially the impact of the plateau is important, as detecting a peak does not, at least in this model, completely fix the dark matter mass.Thus a possible neutrino signal from dark matter may differ from usual assumptions.

Conclusion
Neutrino detection experiments have gained a significant increase in sensitivity over the past few years.It becomes therefore interesting to study their ability to try and detect neutrinos originating from DM annihilations in the present-day universe.In this paper, we have examined one model that allows for such annihilations, the BLSSMIS.A particular difficulty so far in probing the phenomenology that follows from this model is that it features DM particles that couple only very moderately to SM particles.Such DM particles typically have a dominant B B L or h B L component.The resulting DM direct detection cross sections for such particles are well below the neutrino-coherent scattering limit, and the accompanying spectra typically involve neutralino-chargino production cross sections at the LHC that are much lower than those found in the MSSM.However, we find that precisely for these spectra, indirect detection by means of neutrino detection experiments can be used to probe them.The scattering of high-mass B B L or h B L -like DM particles in our present-day universe can create ν h ν h (ν h → ν l + X , where X is any other SM particle except one of the lightest neutrinos, ν l ) and ν h ν l final states.We have shown that a typical feature of such annihilation channels is that they predict a plateau region and a single peak in the neutrino-energy spectrum, which are distinct features that can be measured using neutrino telescope experiments.The extent of the peak region and the location of the energy peak are both completely specified by the mass of the heavy neutrino and that of the DM particle.This shows that measurements of the energy spectrum of cosmic neutrinos can provide a clear and unique way of discovering DM.

χ 0 1 dependence 1 ≈ 62 . 5
due to the chargino-mediated annihilation channel.No W or h solutions are present below 103.5 GeV due to the LEP chargino limits,3 thus the lightest chargino will have a mass similar to the lightest neutralino.The B B L points, and to a lesser extent the B points, show a clear higgs funnel, as can be seen in figure1at M χ 0 GeV.Notably, there is no Z-boson funnel, as we tuned the particle-filter algorithm to search for B B L and h B L LSPs, which favour higgs and Z ′ funnels.These funnels arise due to a higgsino component in both the B and B B L LSPs.The B B L neutralino couples to the SM-like higgs via

Figure 1 :
Figure 1: The relic density of the lightest neutralino as a function of its mass.The dominant component of the neutralino is colour coded as B (blue), W (red), h (yellow), B B L (orange), and h B L (purple).The observed relic density of the Planck collaboration [77] is shaded in blue including an uncertainty band of 0.03 to include computational/theoretical uncertainties.The model points shown here pass all constraints of table 2.

Figure 2 : 1 > 1
Figure 2: The size of the B B L (left) and h B L (right) components for the LSPs in all models that have passed the cuts of table 2, in addition to requiring Ωh 2 < 0.15.Colour coding is the same as for figure 1.

Figure 3 :
Figure 3: The spin-dependent (top) and spin-independent (bottom) DM cross sections for a target proton (left) or neutron (right) of all non-excluded neutralino LSP solutions as a function of the LSP mass.The cross section is given in units of picobarn (pb) and the mass is given in units of GeV.The tentative limits from LZ are shown in blue in addition to the projected limits from PICO, Xenon-nT, and DARWIN.Furthermore, the neutrino-coherent scattering floor is shown in shaded blue.The colour coding for the dominant contribution of the neutralino LSP is the same as in figure 1.

Figure 5 :
Figure5: The DM velocity-weighted annihilation cross section (scaled with the square of their relic density divided by 0.12 and their branching ratio) in cm 3 s −1 as a function of the mass of the LSP (in GeV).The upper panels show the dominant component (upper-left) and dominant annihilation channel (upper-right).For the dominant annihilation channel, we use qq to indicate any combination of quarks produced in a DM annihilation process, ll any combination of leptons, X Y any combination of two bosons such that charge is conserved, ν h ν h any combination of two heavy neutrinos, and ν h ν l any combination of a heavy neutrino and light neutrino.The bottom panels isolate explicitly the LSPs that annihilate mostly into ν h ν h (lowerright) and ν h ν l (lower-left).The Fermi-LAT b b gamma-ray limit and ANTARES ν µ ν µ neutrino limits are shown in the upper plots.The bottom plots show the IceCube and ANTARES ν µ ν µ neutrino limits.Furthermore, the projected KM3NeT νν limits are shown in the bottom plots.The cuts of table 2 have been applied to all shown points, and in addition we require Ωh 2 < 0.15.

Figure 6 :
Figure 6: The neutrino spectrum of two main different annihilation channels:χ 0 1 χ 0 1 → ν h ν h (left) and χ 0 1 χ 0 1 → ν h ν l (right).The three different neutrinos ν e , ν µ , ν τ have been colour coded.All spectra have been computed by MadGraph using 100k events.The predicted peak and plateau regions have been shaded.The width of the shade for the peak region has been widened for visibility.