Sensitivity of direct detection experiments to neutrino dark radiation from dark matter decay and a modified neutrino-floor

In this work we analyse the ultimate sensitivity of dark matter direct detection experiments to dark radiation in form of SM or semi-sterile neutrinos. This flux-component is assumed to be produced from dark matter decay. Since dark radiation may mimic dark matter signals, we perform our analysis based on likelihood statistics that allows to test the distinguishability between signals and backgrounds. Given the previous bounds from neutrino experiments, we find that xenon-based dark matter searches will not be able to probe new regions of the dark matter progenitor mass and lifetime parameter space when the decay products are SM neutrinos. In turn, if the decay instead happens to a fourth neutrino species with enhanced interactions to baryons, DR can either constitute the dominant background or a discoverable signal in direct detection experiments. In the former case, this lifts the “neutrino floor” for xenon-based experiments.


Introduction
The neutral current-induced coherent neutrino-nucleus scattering process [1,2] once inspired the conception of dark matter (DM) direct detection experiments [3][4][5]. Neutrinos with energies up to several hundred MeV elastically scatter on atomic nuclei with a cross section that is approximately enhanced by the squared number of neutrons, N 2 . Similarly, the spin-independent scattering of weakly interacting massive particles (WIMPs) is enhanced by the square of the atomic number, A 2 , boosting the prospects of observing an atomic recoil signal from DM in ultra-low background detectors with keV energy thresholds. Whereas DM has yet to be directly observed in the laboratory, the very process of coherent neutrino-nucleus scattering that once started the a e-mail: marco.nikolic@oeaw.ac.at (corresponding author) b e-mail: suchita.kulkarni@uni-graz.at c e-mail: josef.pradler@oeaw.ac.at field, may also be the defining process in closing the window of opportunity in our direct searches for electroweak-scale DM. Neutrinos produced in the sun, in the atmosphere or in supernova explosions, among other sources, constitute a steady flux that cannot be shielded and, given enough observation time, detector volume, and detection sensitivity will eventually be seen as an irreducible background in DM direct detection experiments. This limits the ultimate sensitivity to discover DM of mass m χ and nucleon cross section σ n , and the combination of parameters where this occurs is conventionally referred to as the "neutrino floor" [6][7][8].
The previous years have seen steady advances in increasing the sensitivity of direct detection experiments. Besides a tremendous effort that is underway and aims at developing and operating ultra-low threshold detectors-see [9] and references therein-the classical WIMP detectors have now gone beyond the ton-year mark in exposure, reaching a sensitivity of σ n 10 −47 cm 2 and better at a DM mass m χ of several tens of GeV. Full exposure results have been reported from liquid xenon experiments LUX [10], PANDAX-II [11], and XENON1T [12], followed by results from current liquid argon detectors DEAP-3600 [13] and DarkSide-50 [14]. The next generation in these experiments, XENONnT [15] and LZ [16] is already under construction and/or commissioning and they will be sensitive enough to see a small number of neutrino background events. Finally, exposures of several hundred ton-year may be achieved with the respective liquid xenon and argon detectors DARWIN [17] and DarkSide-20k [18], and their reach in (m χ , σ n ) will be limited by the neutrino floor.
More than 40 years after its prediction, coherent neutrinonucleus scattering has finally been observed by the COHER-ENT collaboration using accelerator-based neutrino beams [19,20]; efforts to detect the process using reactors are underway [21,22]. These measurements provide valuable new insight into the interactions of Standard Model (SM) neutrinos with the constituents of atomic nuclei, thereby constrain-ing non-standard interactions to quarks and the presence of new forces. In the context of DM direct detection, the most important neutrino source is the sun and beyond-SM neutrino physics utilizing these fluxes has been explored in [23,24]; for more recent works see [25][26][27][28][29][30][31][32][33][34] and the review [35]. A common theme in many of these studies is that new interactions of neutrinos will modify, and typically elevate the standard neutrino floor by within a factor of a few, when the new physics is subjected to complementary constraints.
In this work, we consider a principal alternative option. Rather than modifying neutrino interactions per se, we shall primarily study the presence of new neutrino fluxes, their detectability and their influence on the DM detectability. Concretely, we consider DM decay as a source of SM neutrinos; only in a second step we shall also consider the possibility of new interactions in an extended neutrino sector. Substantial fluxes of these neutrinos will originate from DM decay within our own galaxy as well as globally, by the cosmological decay of DM. Indeed, it is entirely possible that the Universe is filled in significant number with relativistic particles, dark radiation (DR), that may have escaped detection to date. Taking the present energy density in dark matter, ρ DM , as a calibration point, DR may contribute as much as several per cent, ρ DR 0.1ρ DM , while still being allowed by gravitational cosmological probes [36]. When compared to the present number (average energy) of cosmic microwave background (CMB) photons, n CMB ( E CMB ), this implies that Hence, at the expense of having much less DR quanta than CMB photons n DR n CMB , their typical energy may be significantly larger, E DR E CMB . If the energy is in the several tens of MeV ballpark, DR neutrinos induce keV-scale nuclear recoils in direct detection experiments that are principally detectable, and, by extension, also alter the standard expectations for the neutrino floor in DM detection.
In fact, the neutrino energy range between 15-100 MeV is of particular interest because it is a window of opportunityframed by solar and atmospheric neutrino fluxes at the respective low-and high-energy ends-to search for the diffuse supernova neutrino background (DSNB) [37,38]. Its non-observation to-date puts a limit on the flux of electron antineutrinos φ(ν e ) < 3 /cm/s [39], and an upper limit on a cosmological ν e flux from DM decay in this window has been established with Super-Kamiokande (SK) data in [40]; see also [41]. 1 In this work, we shall consider that (a component of) DM decays into neutrinos ν but not anti-neutrinosν. This possibility has been analyzed in some generality in a previous work by some of us [45], where constrains on the combination of lifetime and progenitor mass from neutrino and direct detection experiments have been derived. 2 In this work we explore in detail the DR detectability and the consequences on the DM neutrino floor in the allowed regions of parameter space.
The paper is organized as follows: in Sect. 2 we establish the principal DR fluxes from DM decay, in Sect. 3 we introduce the neutrino-induced event rates at direct detection experiments and list the standard neutrino fluxes, in Sect. 4 we establish the statistical tools for quantifying the principal reach for DM or DR discovery. The main results are then presented in Sect. 5 before concluding in Sect. 7.

Dark radiation from DM decay
We consider the possibility that non-thermal DR is made of neutrinos, either from the Standard Model (SM) or of a new type, that originate from the decay of an unstable DM progenitor X . Due to cosmological constraints [36,47,48], we restrict ourselves to DR scenarios, where only a fraction f X of the total DM abundance injects monochromatic neutrinos via two body decays of a lepton-number carrying Majorontype scalar relic X → νν; a description of such asymmetric model is provided in the original paper [45]. As we are principally interested in neutral current processes, the flavor evolution of injected neutrinos is of no relevance.
The DR flux arriving at the Earth is by and large a combination of two components, the galactic flux ν,gal and the extra-galactic flux ν,e.gal . Assuming no directional sensitivity, for a DM particle X with lifetime τ X and mass m X decaying within the Milky Way, the differential galactic flux is given by, where t 0 = 13.787 ± 0.020 Gyr is the age of the universe, N ν = 2 is the number of neutrinos in final state, E in = m X /2 is the injection energy, r = 8.33 kpc is the distance between the observer at the Earth and the galactic center, ρ = 0.3 GeV/cm 3 is the local DM density and J dec ≈ 2.19 is the angular averaged J-factor obtained from an NFW profile [49]. Compared to t 0 , Eq. (2) probes the amount of DM that is decaying "today." For staying in line with [45], we largely present results for f X = 0.1, but explore smaller, cosmologically better compatible fractions as well.
In contrast, DR that arrives from cosmological distances probes the decaying DM fraction at a time t dec ≤ t 0 , or, in terms of redshift, at z ≥ 0. Hence, the extra-galactic flux is assembled by contributions from all redshifts. To estimate this flux, remember that neutrinos arriving at the Earth with relativistic energy E ν = | p ν | were emitted with higher energy E em ν = E ν (1 + z). The energy-differential extragalactic flux is then obtained via a redshift integral, which for monochromatic injection can be resolved [45], where α = E in /E ν ≥ 1. We set the density parameters dm h 2 = 0.12, M = 0.315, = 1 − M consistent with the CDM model of a flat Universe. H 0 = 100h km/s/Mpc and ρ crit = 3H 2 /(8π G) are the Hubble parameter and critical density at the present time with h = 0.674 [50]. 3 Considering DM lifetimes such that its decay proceeds after matter-radiation equality, we are allowed to neglect the radiation content of the Universe and may find an analytic expression for the cosmic time t (z), where H (z) = H 0 (1 + z) 3 M + is the Hubble rate at redshift z. In this way, a connection between the elapsed time and the redshift is established. Equations (2) and (3) are also applicable for any two body massive final states with replacement E ν → m 2 ν + | p ν | 2 , where m ν , p ν are the mass and momenta of the massive species [45].
The total DR flux X 2ν = ν,gal. + ν,e.gal , a sum of the galactic and extra-galactic component, is obtained by integration over E ν . The galactic flux reads, One may in fact also obtain a useful analytic expression for the extra-galactic flux for the two-body decay with relativistic 3 There is a well-known current discrepancy between the CMB-inferred value of H 0 and the distance ladder estimates from SNIa [51] (among other, more recent local measurements). Our results only depend mildly on the adopted value of H 0 , but we note that decaying DM scenarios such as the one considered here can alleviate this tension [52,53].
final states (see Appendix A), where ρ crit is the Universe's critical energy density today. In the last relation we show the asymptotic scaling for short and long lifetimes, respectively. For τ X t 0 the coefficient is close to unity within 5% because of a numerical cancellation of the value of a logarithm with its prefactor once M and are plugged in. The extra-galactic flux hence asymptotically approaches a maximum for lifetimes τ X t 0 , and can be parametrised as max ν, e.gal In Fig. 1, we exemplify the DR fluxes originating from galactic (dashed lines) and extra-galactic (solid lines) components for various progenitor masses m X assuming a 10% mass-fraction of decaying DM, i.e. f X = 0.1. As both fluxes are inversely proportional to progenitor mass, for a fixed progenitor lifetime the flux decreases as progenitor mass increases. For large lifetime, both galactic and extra-galactic fluxes fall as 1/τ X , for small lifetime, the galactic flux diminishes exponentially, whereas the extra-galactic flux asymptotes to a constant value N ν f X DM ρ crit /m X . This is manifest from (6) and is easily understood: a complete decay of X at some early redshift z corresponds to a mere transferal from one particle species (X ) to another (2ν) with the comoving number of their sum being conserved. "Early DR" originates from a cosmological DM density that is larger by a factor (1 + z) 3 , hence compensating for the flux dilution by a factor (1+z) −3 from area expansion and increase in the timeinterval of subsequent particle arrivals. For τ X t 0 , the total extragalactic flux becomes hence independent of progenitor lifetime. Nevertheless, given the finite energy thresholds of any experiment, the number of "active" neutrinos in a detector diminishes for ν originating at high redshift, because of their redshifting in energy. The maximum flux (as the sum of galactic and extragalactic contributions) is attained when τ X ∼ t 0 . Finally, it is interesting to note that the galactic flux is a double-valued function in m X : for a fixed progenitor mass, the same flux is attained for two different progenitor lifetimes. As we will see later this feature leads to interesting consequences for the experimental sensitivity to DR. The integrated neutrino flux originating from galactic (dashed lines) and extra-galactic (solid lines) components for different DM masses m X as a function of DM lifetime. The flux has the general 1/m X scaling, is inversely proportional to progenitor lifetime for τ X t 0 and the extragalactic component asymptotes to a constant in the converse limit. The maximum DR flux as the sum of solid and dashed lines is obtained for τ ∼ t 0

Signatures of DR at direct detection experiments
Dark radiation in the form of neutrinos produced by unstable progenitors as described in the previous section can potentially be detected at direct detection experiments at the Earth via neutrino-nucleus coherent scattering. Even in absence of such DR neutrinos, the neutrino-nucleus scattering at direct detection takes place due to neutrino fluxes arising from standard ambient neutrino sources. These include solar, atmospheric and supernova-generated neutrinos. We shall refer to those fluxes as the "standard" ones.
The standard neutrino fluxes adopted in this work together with their standard errors are listed in Table 1. For solar neutrinos we use the ones based on the chemical composition determination of [54]. 4 While solar neutrinos dominate for E ν 20 MeV, atmospheric and supernova neutrinos are otherwise the most important fluxes. There are no measurements of the atmospheric flux below 100 MeV and we use the results from a FLUKA simulation [57]. The DSNB neutrinos originate from Type II supernovae and their largest emission in 4 The results of this work are only mildly dependent on the solar opacity problem: the more recent determination [55] primarily affects the O, N, and F fluxes in addition to a 20% downward shift of the 8 B flux. We note in passing that for the last flux, being the most relevant in this context, the measurement [56] lies in between both low and metallicity determinations of [55] and [54], respectively. all flavours takes place during the Kelvin-Helmholtz cooling phase. The prediction depends crucially on the star formation rate [58] as a function of redshift and in this work we use the analytical fit provided in [59]. The emission spectrum is expected to be thermal [60], with each neutrino component being assigned a specific temperature, T ν e ≈ 4 MeV, Tν e ≈ 5 MeV, T x ≈ 8 MeV (x = ν μ ,ν μ , ν τ ,ν τ ). Finally, we note in passing that reactor-and geo-neutrinos are location dependent and relatively small, and we neglect them in this work.
The second relevant aspect is the differential recoil cross section of neutrinos on atomic nuclei. Within SM, the process is mediated by neutral current interactions for spinindependent scattering and can be written as where G F = 1.1663787(6) × 10 −5 GeV −2 is the Fermi constant, Q W = (4 sin 2 θ W − 1)Z + N is the weak charge of the nucleus and E ν is the neutrino energy. The differential cross section is primarily a function of number of Neutrons (N ) because of a near cancellation in the charge (Z ) dependent part of Q W ; sin 2 θ W ≈ 0.23 is the weak angle. The degree of coherence is given by the Helm form factor F(| q|) [64] where q is the three-momentum transfer to the nucleus. As a second possibility we shall consider the case that X decays into a pair of new neutrinos ν B which interact with baryon number through a new vector particle of mass m V and gauge coupling g B [65]. The effective description amounts to a four-fermion interaction In the following we shall only be concerned with relativistic states and we may take the mass of the new neutrino to zero, m ν → 0, so that | p ν | = E ν . The nuclear recoil cross section is coherently enhanced with atomic number A 2 and then reads [45,65], For the sake of presentation, we shall only consider the case of a mediator that is heavy compared to the typical momentum transfer | q| = √ 2E R m N , implying m V 100 MeV in practice. The phenomenology of this model has been explored in [66][67][68][69][70], and in this work we shall consider values G B G F as a possibility to boost the direct detection phenomenology Table 1 Adopted neutrino fluxes in this work; we follow [61][62][63] in their compilation of the fluxes (original works are referenced in the main text). Beside the fluxes, the endpoint or the maximally used energy together with the corresponding maximal nuclear recoil energies on a xenon target are given of DR. 5 Current limits on G B arise from rare meson and Zdecays such as K → π V , B → K V and Z → V γ [68,71], respectively, monojet searches [72,73], as well as from neutrino facilities [67,74]. For example, for m V = 300 MeV the most stringent guaranteed constraints come from Mini-BooNe, requiring G B 500G F . The limits significantly tighten by an enhancement of longitudinal V emission when baryon-number is broken by the chiral anomaly, essentially ruling out G B > G F [68,71]. Although those limits are fairly robust, they may also be evaded depending on the anomalycanceling UV completion of the model [68], and, G B > G F remains a possibility that we entertain as the most optimistic case in this work. The differential recoil rate introduced by neutrinos of source i is therefore given by 5 Besides DM decay, ν b can be generated from oscillation or through direct production. The former contribution is rendered unimportant is the new (small) mixing angle or by virtue of an astronomical oscillation length that prevents their appearance in terrestrial or solar fluxes [24]. The direct production of ν b may likewise be neglected for the purpose of this study: any atmospheric ν b -yield competes with the QCD production of mesons and is much smaller; any thermal solar ν b flux (e.g. through proton-proton bremsstrahlung) is undetectable on the account of keV-scale energies, as is its SM counterpart. Finally, a DSNB in ν b is an interesting possibility, but the final spectrum is expected to be suppressed both in energy and overall flux because of a very efficient trapping inside the supernova. Out of these reasons we expect little to none interference with the DR fluxes and therefore neglect those additional contributions.
where E ν,min = √ E R m N /2 is the minimum neutrino energy to produce a recoil on a target of mass m N ; E max ν is the maximum energy neutrinos of source i can have. In case of DR, E max ν is given by endpoint energy of the source E in = m X /2. The number of target nuclei per unit mass of detector material is denoted by N T . The total rate in a detector will then be the sum over all isotopic compositions of relevant elements and over all neutrino sources i. We shall denote by d R ν (E R )/d E R the total rate over all "standard" sources and by d R X 2ν (E R )/d E R the DR-induced non-standard recoil rate. The latter carries two contributions, from galactic and extragalactic fluxes, respectively. An analytic expression for the differential recoil rate for the extragalactic flux can be obtained and is given in Appendix A.
The maximum nuclear recoil energy from DR is derived from kinematics and depends on the mass of the progenitor, where E max ν = E in is the neutrino energy at injection. It is worth noticing that the end point of recoil energy is directly proportional to the square of progenitor mass. The expected recoil rate is then given by the integral over the recoil energy, where E R thr is the threshold of the detector.
As an example, consider the recoil rate induced by a progenitor of 100 MeV mass decaying to SM neutrinos in a liquid xenon detector with negligible nominal 1 eV threshold. Notwithstanding a general dependence on recoil energy, saturating the possible DR fluxes, we find that the differential rate is bounded from above by ) .
This implies for the total detectable rate, and demonstrates that for a given progenitor mass, there is a natural ceiling on the event rate when considering DR in form of SM neutrinos. In Fig. 2, we show exemplary differential recoil rates in xenon introduced by DR in form of SM neutrinos for a progenitor mass of 60 MeV in comparison with the differential recoil rate due to standard neutrino fluxes. The progenitor lifetime takes on values from 10 −3 Gyr to 10 4 Gyr. Heavier progenitor particles introduce larger recoil energies at direct detection experiments, as the recoil endpoint energy is controlled by m X . Generally speaking, m X 30 MeV is required to exceed solar neutrino-induced events-in particular the ones induced by 8 B neutrinos-because of the principal limitations in the magnitude of incident DR flux.
The latter point is further exacerbated by the fact that for a percent-level decaying fraction of DM we find that there is no combination of (τ X , m X ) such that the total DR-induced event rate in xenon is larger than the one from 8 B neutrinos, While it is well known that approximately 6 GeV mass WIMP recoils mimic 8 B neutrino recoils, we find that DR induced neutrino recoils from 2-body DM decay can neither mimic nor exceed the latter. This puts a principal limitation on the detectability/influence of DR in direct detection experiments: any set of prospective parameters must be such that the DR flux has a component that exceeds the 8 B flux in energy, hence pointing to m X 30 MeV. 6 Before discussing the employed statistical framework, for completeness, we close this section by reviewing the DM (χ ) signal hypothesis, which has not been discussed so far. The total event rate of DM recoils at the direct detection experiments, integrated over the differential recoil spectrum In the second equality, we make the assumption of a standard spin-independent DM-nucleus interaction of contact type. The local DM mass density is ρ χ 0.3 GeV/cm 3 and m χ (μ χ,n ) is the DM (reduced DM-nucleon) mass and σ n is the DM-nucleon scattering cross-section with equal coupling to protons and neutrons. At last, f ( v ) is the normalized galactic Maxwell-Boltzmann velocity distribution with v 0 = 220 km/s, truncated at the escape velocity v esc = 544 km/s and boosted into the detector frame with v lab = 232 km/s. Detector specifics are only entering via the recoil threshold energy E thr and assume 100% efficiency of detection for E R ≥ E thr , ignoring finite energy resolution effects. Furthermore, for our purposes it is sufficient to integrate the 6 More precisely, if the statistical error bar on the 8 B-induced events becomes smaller than the DR-induced event rate, one may actually subtract the 8 B background and DR becomes detectable. This is seen by the extra island at m X 10 MeV in left Fig. 5 for the futuristic exposure of 100 ton-year (see below). differential recoil rate up to a recoil energy of 100 keV. It should be noted that μ χ , μ ν , and μ X 2ν represent recoil rates per detector mass and live-time and need to be multiplied by detector exposure ε to get the total number of observed events.
We now introduce the formalism to make quantitative statements on sensitivity of direct detection experiments to DR and its influence on the observation of a potential WIMP signal.

The principle reach of direct detection experiments
To forecast sensitivity of an experiment, there are two seemingly similar, yet distinct questions to address: the first regards the ability to exclude the presence of a signal and the second regards the ability to discover the signal. The exclusion and discovery exercises can in turn be performed either assuming no backgrounds or in the presence thereof.
The purpose of this work is to see to what extent direct detection experiments can discover DR. Furthermore, we like to quantify the reach of direct detection experiments to a DM signal in presence of DR in addition to the standard neutrino background. We therefore deal with two different signal and background hypotheses, one which involves DR recoils as signal and standard neutrino sources as background and one which involves a DM signal and DR plus standard neutrino sources as backgrounds. In each case, we are interested in understanding the exclusion potential and discovery reach. The methodology we outline below closely follows previous work on the neutrino floor in direct detection [7]; for an alternative approach see e.g. [75].

DR exclusion potential of direct detection experiments
In this part, we set up the formalism to obtain the potential of direct detection experiments to exclude a DR signal for as long as backgrounds remain absent. This situation can be dealt with easily, and we outline it below; the signal hypothesis in the presence of background will be dealt with using the likelihood approach as in Sect. 4.2.
For a single channel counting experiment, the probability to observe n events, when λ events are expected is given by the Poisson distribution Pois(n|λ) = λ n e −λ /n!. In general, λ can be the expected number of signal events εμ X 2ν , background events, or their sum. For zero observed events, n = 0, and zero expected background, the 90% CL upper limit on λ is obtained from a p-value p λ = 0.1, corresponding to a 10% probability that the outcome of an experiment is at least as extreme as observed. It is then related to the confidence level via p λ ≤ 1−CL. Solving p λ = Pois(0|λ) for λ yields λ = 2.3 as usual.
In a background-free experiment with threshold energy E thr and exposure ε, the 90% CL exclusion contour in the (m X , τ X )-plane is found by evaluating (12) under the condition εμ X 2ν = 2.3. The strongest exclusion that can be obtained in such background-free scenario, is the one where the exposure grows to the level that the first irreducible neutrino background event is seen, Here we have highlighted the threshold-dependence of the required exposure, ε(E thr ); the neutrino-induced recoil rate is given by the sum over all "standard sources" in Eq. (10). Later we shall also consider the possibility that DR is present as additional "background," and we wish to understand how it impacts the standard literature results on the DM neutrino floor [7]. The condition for when background starts limiting the DM in the presence of DR search reads in an obvious modification to (15).

DR discovery potential of direct detection experiments
While the procedure described above is applicable in the zero background scenario, the upcoming/ongoing ton scale direct detection experiments with large exposures will not remain background free. These backgrounds which originate from standard neutrino recoils suffer from uncertainties in the measured neutrino fluxes (see Table 1) and which need to be taken into account when charting out the discovery reach of an experiment.
We assume that observed spectra of events are composed from signal and background sources. As such, they enter the generalized Poisson likelihood, where the sum over sources α depends on the hypothesis under question; model parameters such as m X and τ X are part of the vector θ . In our simulations, the individual recoil events E R i are random numbers, a total N α of them drawn for each source α from the probability density function (PDF) μ −1 α d R α /d E R . Note that N α itself is a Poisson random number with mean value μ α . Finally, N obs = α N α is the total number of observed events.
When a source of events is declared a background under the considered hypothesis their associated fluxes become a nuisance parameter. Concretely, we model the flux uncertainties as Gaussians with mean φ β and variance σ 2 β , Eur. Phys. J. C (2022) 82:650 Here, the product is over all background sources β and the random variables φ β are part of θ. For the standard neutrino sources the mean values φ β together with their standard errors are listed in Table 1. When DR is considered as a background, we take φ X 2ν = ν,e.gal + ν,gal and assume σ X 2ν = 0.3. The assignment of a 30% uncertainty is somewhat arbitrary, but meant to be conservatively reflective of overall astrophysical uncertainties concerning DM-induced DR fluxes. Taken together, the total likelihood becomes, For a given threshold energy E thr and exposure ε and for every point θ in the parameter space, we then generate several hundred mock realizations to build a statistical sample. Each of these realisations we refer to as an "experiment". We now spell out the signal and background hypotheses that we are considering in this work: (1) DR signal with standard neutrino backgrounds: The first step is to explore the DR discovery potential of direct detection experiments. Hence, we treat the DR recoils as signal in the presence of standard neutrino backgrounds but in absence of DM-induced events. The null hypothesis is H 0 : σ n = 0 and τ X → ∞ and the alternative hypothesis becomes H 1 : σ n = 0 and finite τ X . We therefore take the sums in Eqs. (16) and (17) α = X 2ν, ν j and β = ν j , respectively. In the baryonic neutrino scenario, we fix G B to an exemplary value.
(2) WIMP signal/no DR background: As a validation procedure we may also consider the complete absence of any DR. This allows us to recover the standard result on the neutrino floor in the (σ n , m χ )-plane [7]. In order to discover a WIMP signal, the background only hypothesis H 0 : σ n = 0 must rejected; the alternative hypothesis is H 1 : σ n > 0. We take α = DM, ν j in Eq. (16) and β = ν j in Eq. (17) to minimize the likelihood ratio (19).
(3) WIMP signal with DR background: The final alternative is to switch on DR as an additional source of background and study the influence on the WIMP discovery region. The background only hypothesis remains H 0 : σ n = 0 and the alternative hypothesis is again H 1 : σ n > 0, but now we take the sums in Eqs. (16) and (17) to include DR, i.e. α = DM, X 2ν, ν j and β = X 2ν, ν j , respectively. We fix progenitor mass and lifetime for this procedure.
Once H 0 and its alternative H 1 are chosen, one may then proceed in standard fashion and use the negative loglikelihood as test statistic for each experiment [76], whereˆ θ maximizes the likelihood under the background-only hypothesis H 0 andˆ θ maximizes L for signal plus background, H 1 . The distribution of q under H 0 asymptotically follows a χ 2 -distribution with one degree of freedom as per Wilk's theorem [77]. The p 0 -value then characterizes the probability that the background-only hypothesis H 0 is excluded if the qvalue is greater than the "observed value" q obs , where f (q| H 0 ) is the PDF of q under H 0 . We obtain the latter by Monte Carlo generation of mock data. In terms of significance Z = √ q, given in units of standard deviations, the p 0 -value is obtained by the cumulative standard normal distribution (x), For a discovery with 3σ significance the corresponding p 0value is given by p 0 = 0.00135. If data is observed in the critical region Z ≥ Z obs the hypothesis H 0 is rejected. For the 3σ significance that we are interested in, Z obs = 3. We have verified, that our generated background data satisfies In turn, the probability β for H 1 being rejected is defined by, where f is now the PDF of Z under H 1 and Z 90 represents the significance that can be at least obtained 90% of the time in an experiment with data generated under H 1 . Conversely, H 1 is accepted at a confidence level 1 − β. Assuming a confidence level of 90% the alternative hypothesis H 1 is excluded for β = 0.1. It is thus possible to find a value of signal with 1 − β = 0.9, which leads to Z 90 ≥ 3. For this case, an experiment has a 90% probability to detect at least a 3σ signal. Concretely, we generate between 250 and 500 Monte-Carlo data sets for the signal plus background hypothesis H 1 for each value of progenitor mass m X or DM mass m χ and for each lifetime τ X or cross section σ n . These data sets determine the significance distribution. The procedure is repeated for a range of lifetimes or cross sections until we obtain Z 90 = 3 for β = 0.1 which corresponds to a 3σ discovery potential at 90% CL. We repeat this procedure until the entire parameter space is mapped out.
The formalism described above is not specific to any particular detector. As was shown above, DR fluxes from DM decay can compete with standard solar neutrino fluxes only above the 8 B endpoint energy. Dark radiation introduces changes to the standard picture in the recoil energy region where the large WIMP detectors have ample sensitivity and ultra-low thresholds are of little benefit. In the following, we will hence explore the formalism on the basis of ton-scale xenon experiments like the ones mentioned in the introduction.

DR exclusion limits
We first establish the exclusion potential for DR in the presence of standard neutrino backgrounds. In the remainder of the paper we shall focus on the example of a liquid xenon detector. The operational procedure is then as follows: given the presence of standard neutrino fluxes, we first compute the DR limit at a fixed threshold by finding the exposure for which one background event as explained in Sect. 4.1 is seen. In a second step, the threshold is varied and the hull of all exclusion curves is obtained. The resulting contour in the (m X , τ X ) plane then represents the optimal DR exclusion potential of the experiment when irreducible backgrounds are not dealt with. Figure 3 shows the results of the procedure, i.e., the exclusion potentials for DR in form of SM neutrinos (left) and baryonic neutrinos (right) for threshold energies ranging from 0.01 keV up to 50 keV. In the region τ X > t 0 they are determined by the galactic and extra-galactic components whereas in the region τ X < t 0 only the extra-galactic flux contributes. The prospective regions are hence centered around t 0 : for τ X t 0 the flux is too small and for τ X t 0 the extra-galactic flux becomes too soft. It is interesting to note that in the right panel of Fig. 3 below the minimum assumed threshold of 1.2 keV, DR in from of SM neutrinos cannot be excluded, as their induced events are superseded by solar neutrinos, especially the ones originating from the 8 B reaction. On the flip side, as the threshold increases, the experiment becomes more sensitive to higher energy recoils generated by larger m X . Hence, in contrast to the WIMP case, and because of the general limitations in the DR flux in form of SM neutrinos, ever lower thresholds as they are sought in many direct detection experiments are not a beneficiary factor.
So far, for simplicity, we assumed f X = 0.1 throughout. At this point, it is worthwhile exploring the impact of smaller values of f X that remain cosmologically unconstrained. This is shown in Fig. 4, and one should bear in mind that the exclusion limits are iso-contour lines of 2.3 DR events at the optimal threshold-dependent exposure where a SM neutrino background starts to be seen. In the large lifetime region τ X t 0 , f X only changes the overall flux normalization: a decrease in f X hence demands a decrease in τ X for a fixed progenitor mass m X to attain the same number of events. In the small lifetime regime τ X t 0 , events are due to the extra-galactic component. Here the position of curves is rather determined by the fraction of flux that has not been redshifted below the relevant threshold value that is associated with the model-point (m X , τ X ). Overall, it is a combination of factors that implies that the discovery reach is more slowly deteriorating with decreasing f X . Finally, in the right panel, one additionally observes the strong decrease in sensitivity for τ X ∼ t 0 with falling f X . It happens because the 8 B background becomes visible before 2.3 DR events can be attained with smaller f X at the associated exposure; this is the situation all along in the SM case (left panel) which does not have such "nose" feature even for f X = 0.1. Similar conclusions with respect to a decrease in f X will also hold for the discovery potentials discussed in the next section.
The exclusion potentials have to be put in perspective with the current sensitivity of multi-ton scale neutrino experiments. This has been explored previously in Ref. [45] and the gray shaded regions in Fig. 3 show the excluded regions in the DR progenitor plane (m X , τ X ) from Super-Kamiokande and Borexino. The constraints included in the left panel of the Fig. 3 were derived from measurements of solar-and atmospheric fluxes [78,79], and from searches of DSNB neutrinos [39]. In the right panel of Fig. 3 the constraint was obtained by re-purposing a Borexino solar-axion search [80] to ν B -induced elastic proton recoils; for further details on both cases see [45]. 7 The inside of the curves labeled according to the assumed threshold correspond to a combination of parameters where we find that DR sourced by a DM progenitor is discoverable in a direct detection experiment.
In summary, except for relatively small regions located in the quadrant m X 100 MeV and τ X > t 0 , neutrino experiments carry stronger current constraining power to SM DR than DM direct detection experiments will every have.
On the other hand, when we are to entertain new physics in the DR itself (in form of ν B ), direct detection sensitive to lower masses m X 20 GeV is gained as the threshold drops below 1 keV. But even with higher threshold values, and expanded region in τ X at a given m X can be covered.

DR discovery potential
We now turn our attention to the prospects of making a discovery of DR in the presence of irreducible backgrounds which are to be faced in the coming generation of direct detection experiments. Charting out the discovery regions will be achieved using the likelihood approach outlined in Sect. 4.2. Similarly to the previous section, we attempt at discovering DR recoils in presence of standard neutrino backgrounds. Importantly, such exercise requires specification of an experiment's exposure and threshold and is hence specific to these assumptions. In the following we choose ε = 1 ton-year and E thr = 0.1 keV to normalize on current collected data sets combined with an aggressive assumption on nuclear recoil threshold. Figure 5 shows the result. We consider exposures up to 100 ton-year for a threshold of 0.1 keV for DR in form of SM neutrinos (left panel) and baryonic neutrinos (right panel). For the extreme case of 100 ton-year exposure, reflective of the sensitivity of future multi-ton xenon detectors such as DARWIN [17], we switch to a binned likelihood approach with N bin = 50. The limiting factors for the reach can be understood by analysing the most dominant background fluxes. In the bulk of the (m X , τ X ) parameter space, hep neutrinos constitute the dominant background to DR. This is because hep neutrinos generate a large range of recoil energies with appreciable flux. Around 20 MeV progenitor mass, the 8 B neutrinos become the limiting factor. As can be seen in the left panel of Fig. 5, it is not possible to make a 3σ discovery of DR in form of SM neutrinos, unless a 100 ton-year exposure is assumed (and paralleling progress with neutrino detectors is neglected.) When ν B are invoked The distribution of maximum likelihood parameters σ n , m χ under the DM hypothesis when they originate from DR with various combinations of (τ X , m X ) in form of SM neutrinos (left panel) and baryonic neutrinos (right panel) instead. The gray points as labeled show excluded combinations of (τ X , m X ) for an exposure of 1 ton year and threshold of E thr = 2.5 keV. The colored points show allowed combinations of (τ X , m X ) for an exposure of 100 ton year and threshold of E thr = 0.1 keV (right panel), they can still be discovered in a wide mass range with only 1 ton-year exposure.

Dark matter vs. dark radiation
Here we shall first consider the scenario in which a signal of DR origin is interpreted as DM. For this purpose, we generate one hundred mock-data sets where we inject a DR signal on top of the usual neutrino background. In order to study how such a false positive DM signal is seen, we fit the mockdata by a binned version of the DM-only likelihood for α = DM and find the best fit parameters σ n , m χ . The results of this procedure are shown in Fig. 6 where in the left (right) panel we consider SM (baryonic) neutrino DR. The gray dots and crosses show otherwise disfavored parameter points in (τ X , m X ) that are within the discovery reach for a threshold of E thr = 2.5 keV but avoid the large 8 B background. In turn, the colored points and crosses are for allowed parameter combinations, but require a low threshold E thr = 0.1 keV for discovery and hence face a large solar neutrino background. The general trend in both panels is, that for larger progenitor masses m X , mock-data is fitted by larger DM mass m χ . The lifetime dependence of DR is primarily imprinted into a linear scaling of the cross section σ n and the highest fitted cross section is given for τ X ∼ t 0 , when the DR flux is maximal. For some combinations (τ X , m X ) for DR in form of baryonic neutrinos, the cross section is even higher due to the enhancement A 2 G 2 B in the neutrino cross section. For SM DR (left panel) the most important feature in the misidentification is that when the DR flux diminishes, the fitted points start tracing out the boundary of the neutrino floor. When the smaller threshold is considered, all fits indeed center in on the 8 B tip, at (σ n , m χ ) ∼ (4 × 10 −45 cm 2 , 6 GeV), of the neutrino floor. This shows that although SM DR can be discovered with 3σ significance in an experiment with exposure ε = 100 ton year as shown in the previous section, it requires a careful statistical analysis that goes beyond the naive fits presented in this section, for which one would rather make the erroneous conclusion of only facing the neutrino floor. For the baryonic case (right panel) the prospects are better in the sense, that a naive fitting would largely point to a new physics signal (even if misinterpreted in term of its origin.) 6.2 DM exclusion potential in presence of DR Next we consider the modifications to the exclusion potential for DM when DR neutrinos become a source of additional background. The procedure is similar to before: given the presence of an anomalous neutrino flux, we first compute the best DM limit at a fixed threshold by finding the exposure for which one background event is seen. Varying the threshold, the minimum value of σ n for each m χ , is obtained.
In the left panel of Fig. 7, the results for various currently allowed (m X , τ X ) combinations are shown. In the right panel, a progenitor mass of 60 MeV, and lifetime of 10 Gyr decaying into ν B with G B = 10G F is shown. The various dashed curves show the best limit on σ n for various assumed thresholds, from 1 eV to 50 keV, as labeled. The blue (red) solid thick line is the minimum of all curves and represents the optimal exclusion potential. The curves are to be compared with the gray line and shaded region which represents the standard exclusion potential.
In the left panel, only modest modifications to the standard result are obtained in the regions between 6−50 GeV WIMP mass. In the right panel, the sizable new coupling boosts the DR-induced nuclear recoil rate above the 8 B neutrino recoil rate, weakening the minimal excludable DM-nucleon cross section by few orders of magnitude once m χ < 50 GeV.
The general trend in both panels of Fig. 7 is that the exclusion potential is shifted to the right off the standard 8 B shoulder. How far that shift goes, depends on kinematics. For a progenitor mass of m X = 60 MeV, DR neutrinos induce recoils that compete with DM only below m χ 50 GeV. The exclusion potential for heavier DM is then only challenged by the standard, more energetic atmospheric fluxes. Only if we were to entertain the possibility of m X 100 MeV into the GeVregime (generally challenged from neutrino experiments), the floor would be lifted by several orders of magnitude in the electroweak scale DM mass regime. GeV, when DR fluxes come to dominate the event rates. The modification of discovery limits are seen up to DM masses of 1 TeV. As can be seen in the left panel, the maximal allowed DR signal is at most comparable to the SM neutrino backgrounds and only mildly alters the standard neutrino floor (gray dashed). In the right panel, because of large statistics, we switch to a binned likelihood approach with N bin = 50. Owing to the increased strength of interaction, the neutrino floor is also raised in the low mass region left of the 8 B shoulder (m χ 6 GeV).
In both panels of Fig. 8, the modification of discovery limits are seen up to the largest shown DM masses of 1 TeV. This is in contrast to the exclusion potentials of the previous section. For the latter, only the number of events for varying threshold (up to 50 keV) enter, but not their detailed spectral shape. In this section, we fix the threshold to E thr = 0.1 keV so that DR spills into the standard induced events at any rate, also affecting the discovery potential for m χ 50, GeV.

Conclusions
In the foreseeable future, DM direct detection experiments will face irreducible background in the form of neutrinonucleus elastic scattering. The responsible neutrino fluxes originate from the sun and to a smaller degree from the atmosphere or from the cosmological SN history. In this work, we first study the detectability of neutrino DR, sourced by a percent fraction of DM that may decay with arbitrary lifetime after recombination. We consider either SM or baryonic neutrinos ν B as DR. The latter designate a semi-sterile fourth species that interacts through gauged baryon number. In a next step, we also explore the impact on the DM discovery and exclusion potential at next generation direct detection searches in the presence of new, anomalous neutrino fluxes.
The DR fluxes principally fall into two categories, either constituting neutrinos that travel cosmological distances or neutrinos that originate from DM decay in the galaxy. Whereas the latter requires τ X t 0 as otherwise the DR-generating sub-component of DM would have already decayed, the former contributes for any τ X and we provide a new closed-form expression for its total flux. In the detector, the neutrinos then interact coherently with the nucleus, with a cross section that approximately scales as (A − Z ) 2 and A 2 for SM and baryonic neutrinos, respectively. The benefit of DR in form of SM neutrinos is that its interactions are known. In turn, DR in form of ν B allows for effective interactions that can be stronger than in the SM, G B > G F , boosting the principal importance of the DR signal.
After establishing the exclusion potential of DR, we then address the more important question of detectability of DR on the example of a ton-scale liquid xenon detector. Here, standard neutrino sources act as the irreducible background. Given that the signal might be salient in the data, we employ a profile-likelihood analysis to reach robust statistical conclusions about the various assumed hypotheses. For various combinations of detection threshold and exposure, we build a statistical sample by Monte-Carlo generating a large set of mock recoil spectra that are subsequently used in the maximization of the likelihood function. Through this numerically expensive procedure we are able to obtain the discovery reach to DR as well as a DR-modified "neutrino floor", i.e. the boundary in the (m χ , σ n ) parameter space for the ability to detect DM with 3σ significance when neutrino backgrounds cannot be rejected.
Concretely, we find that SM DR cannot be discovered in direct detection experiments unless a futuristic exposure of 100 ton-year is assumed. The coverage in new parameter space is modest (and neglects the paralleling advances in neutrino detectors). In contrast, direct detection experiments bear better potential to discover non-standard neutrino interactions of DR, as exemplified in the ν B case, see Fig. 5. Here progenitor masses down to m X = 10 MeV and lifetimes as large as 10 5 t 0 can be discovered.
In a next step, we turn to the potential DR-induced interference for the DM searches. Here we study the confusion potential between DR and DM and establish the DM exclusion and discovery potential in the presence DR in future WIMP searches. For DR in form of SM neutrinos, we find that the best sensitivity to the DM-nucleon cross section is weakened by less than an order of magnitude once the scenario is subjected to existing constraints from SK data. Considering ν B as DR with G B > G F we find that the entire region m χ 50, GeV can be strongly affected, and the DM exclusion or discovery potential can be modified substantially.
With the advance of multi-ton scale direct detection experiments, we are entering the "end-game" of WIMP detection. Neutrinos from standard sources will become an irreducible background and ultimately limit our ability to push for ever smaller event rates and better DM sensitivity. It is hence only timely to ask, what other signals and irreducible backgrounds there could be. Here we investigated the perfect possibility that our Universe is filled with MeV-scale DR that traces back to the instability of (a component of) DM. We find that standard expectations can be altered, and-in the presence of new particles and interactions-in a significant manner. In the latter case, DM direct detection experiments may then be turned into DR detectors instead.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The data can be easily generated by random Monte Carlo experiments as described in Section 4 and consist of random recoil energies from the signal and background hypotheses].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, pro-vide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.

Appendix A: Analytic form of d R X2ν /d E R
For 2-body decays into massless neutrino DR, the energydifferential recoil rate given in Eq. (10) can be integrated analytically. For the extragalactic signal with the flux given in (3) we split the integral over neutrino energy into the two terms that comprise the cross section (8), where E in = m X 2 is the injection energy and E min = E R m N 2 is the minimum energy to produce a recoil E R . The integrals are solved by substitution, and we obtain the expression, where The total flux X 2ν in Eq. (6) is obtained from the first two terms in the bracket of (25) by setting x min = 1. For the galactic component the integral with the flux given in (2) can be evaluated directly, The DR rate is then the sum of both, galactic and extragalactic contributions. The baryonic neutrino rate is calculated in the same way replacing Q 2 W G 2 F /4π by A 2 G 2 B /2π .

Appendix B: Details on discovery and exclusion potentials
Here we show some examples of exclusion and discovery limits, where SM neutrinos may mimic potential DR and DM signals in direct detection experiments. The reach of an experiment for a signal depends on the shape of the recoil rate and the number of events in relation to backgrounds. In this appendix we study cases where DR neutrinos exhibit some degree of similarity with DM-induced events.
In the left panel of Fig. 9 we show two examples where DM and DR yield similar differential recoil spectra for two different progenitor masses 50 MeV and 1 GeV and DM masses as labeled such that the recoil endpoint energies match. For the lighter progenitor-and which is most prospective with respect to existing constraints (see above)-the DM and DR spectra are well distinguishable, with DR-induced events exhibiting a harder spectrum. For GeV-scale progenitors, the spectrum is cut off by the nuclear form factor |F(q)| 2 explaining thee similarity in the rate with an electroweakscale mass DM particle. We note in passing, that in the latter case, additional inelastic nuclear channels are possible as the CM energy is in the tens of MeV regime.
In the right panel of Fig. 9 we follow our usual procedure and simulate several thousand Monte-Carlo experiments for each neutrino source ν j and calculate the distribution of best parameters (m X , τ X ) by minimizing the likelihood from Eq. (16) for α = X 2ν. The discovery limit to SM DR neutrinos for an exposure of 1 ton-year and threshold E R thr = 0.1 keV is shown. In addition, the dots as labeled are indicative of the standard neutrino fluxes that limit a further extended reach.