Signatures of Dark Radiation in Neutrino and Dark Matter Detectors

We consider the generic possibility that the Universe's energy budget includes some form of relativistic or semi-relativistic dark radiation (DR) with non-gravitational interactions with Standard Model (SM) particles. Such dark radiation may consist of SM singlets or a non-thermal, energetic component of neutrinos. If such DR is created at a relatively recent epoch, it can carry sufficient energy to leave a detectable imprint in experiments designed to search for very weakly interacting particles: dark matter and underground neutrino experiments. We analyze this possibility in some generality, assuming that the interactive dark radiation is sourced by late decays of an unstable particle, potentially a component of dark matter, and considering a variety of possible interactions between the dark radiation and SM particles. Concentrating on the sub-GeV energy region, we derive constraints on different forms of DR using the results of the most sensitive neutrino and dark matter direct detection experiments. In particular, for interacting dark radiation carrying a typical momentum of $\sim30$~MeV$/c$, both types of experiments provide competitive constraints. This study also demonstrates that non-standard sources of neutrino emission (e.g. via dark matter decay) are capable of creating a"neutrino floor"for dark matter direct detection that is closer to current bounds than is expected from standard neutrino sources.


INTRODUCTION
The dominance of dark matter (DM) and dark energy (DE) in the total energy balance of the Universe is a widely acknowledged nonetheless astounding fact. Precision studies of the cosmic microwave background (CMB) [1] provide evidence that DM was "in place" before hydrogen recombination, while the effects of DE manifest themselves much later, starting from redshift z ∼ O(1) [2]. Given a rather elaborate structure of the visible sector, the Standard Model of particles and fields (SM), it would be a logical imperative to consider somewhat more complicated models of a dark sector, beyond a single particle contributing to DM, and a cosmological constant sourcing DE. Indeed, in recent years there has been some significant progress in studying models of dark sectors [3], which include the possibility of new "dark forces", new massless states (also known as "dark radiation", or DR), and/or relativistic massive dark states that can be produced through late time processes (also known as "boosted dark matter") [4][5][6][7], co-existing with massive DM particles. Both new massless states and massive boosted states in the dark sector can be classified as "dark radiation" in general terms. The presence of such sectors significantly broadens the phenomenology of DM [5,6,[8][9][10][11][12][13][14][15][16][17][18][19][20][21] motivating, in turn, a wider scope for the experimental efforts dedicated to the searches of DM.
An enormous progress in observational cosmology has resulted in a very sensitive constraint on the number of extra degrees of freedoms that remained radiation-like during the CMB epoch. The Planck collaboration has reported a stringent constraint, that phrased in terms of neutrino-like radiation species reads as [1]: where ρ DR is the energy density in additional dark radiation. With further refinements [22,23], one can show that the actual 2σ limit on the deviation from the SM prediction for N eff is ∆N eff ≤ 0. 39. This bound has a wide degree of applicability, and is most effectively used to constrain models with "early" DR, or models with extra light degrees of freedom that were in thermal contact with the SM, but decoupled at some point in the history of the early Universe. In that case there are also extra bounds provided by big bang nucleosynthesis (BBN), that can be often cast in terms of the same parameter N eff [24][25][26]. However, constraint (1) would not be applicable to models where DR is created much later than the CMB epoch. For example, recent decays of a sizable fraction of DM into dark radiation are allowed, and, moreover, ρ DR can be much larger than ρ γ today.
In this paper, we are interested in the late generation of ρ DR with the following properties: the number density of DR particles is smaller than that of CMB photons, while the kinetic energy on average is much larger than the energy of individual CMB quanta, In the last relation, we require that the amount of dark radiation does not exceed 10% of the dark matter energy density, in accordance with recently updated constraints [27]. Such set of inequalities leaves, of course, a lot of freedom for what DR can be, but restricts a number of possibilities for how the non-thermal DR can be created.
In this paper, we will consider a scenario where the dark sector mediates some DR-SM interaction to be specified below. The new interactions allow to probe DR directly via its interaction with nuclei and electrons rather than via its contribution to the Universe's energy bal-ance (through ρ DR .) Cosmic SM neutrinos are, in some sense, the example of interacting dark radiation. Remnants from the Big Bang, they have a very small energy today E ν ∼ m ν , and their direct detection via weak interactions represents a huge experimental challenge [28]. (There is, of course, plenty of evidence for cosmic neutrino weak interactions in the outcome of BBN.) It is also known that neutrinos have other cosmic components, including the one generated by global activity of supernova (SN) remnants [29]. The search for the diffuse SN neutrino flux is a challenging endeavor, driving some developments in solar neutrino detectors [30].
Neutrinos provide a small-but in future importantbackground for the searches of weakly interacting massive particles (WIMPs) in direct detection experiments [31][32][33]. 1 The solar neutrino flux is the largest component for such a background, but given a relatively low energy cutoff to its spectrum, the generated nuclear recoil is rather soft. Above a neutrino energy of 18 MeV only diffuse supernova and atmospheric neutrinos constitute the primary components of the neutrino spectrum. It has to be emphasized that there is no direct observation of the diffuse supernova flux yet, and only an upper limit exists on the flux of electron antineutrinos in the 15-30 MeV window [37], Φν e < 3/cm 2 /s. Above 30 MeV, the atmospheric neutrinos start to be dominant, with measurements above 150 MeV energy [38], and a total flux in the ballpark of a few/cm 2 /s. The above window is quite important because it corresponds to a momentum transfer scale that is associated with the optimum sensitivity region of dark matter detectors, such as large xenon-based detectors [39][40][41]. Therefore, if there is an additional neutrino or neutrino-like component of dark radiation, the bottom of the direct detection "neutrino floor" can be closer than expected. 2 Given the huge amount of efforts devoted to the scaling-up of the WIMP direct detection experiments, one should investigate possible signals from additional hypothetical components of the neutrino flux, and from DR in general. Apart from creating a new signal competing with WIMPs, the search for DR is interesting in its own right, and can be done in parallel to WIMP searches. It is well known that direct detection experiments often provide the most sensitive tool for broad classes of physics with 1 eV to a few 100 keV energy release (see, e.g., Refs. [45][46][47][48][49][50][51] for a sample of rep- 1 The idea to study the coherent scattering of neutrinos on nuclei [34] pre-dates the DM direct detection idea [35], and only recently has it been observed with neutrinos sourced by meson decays [36]. 2 An upper limit on a cosmological flux of SM neutrinos from DM decay in the ∼15-100 MeV mass window using Super-Kamiokande data has previously been established in [42]; see also recent Ref. [43]. The ensuing lower limit on the DM lifetime is driven by theνe flux-component; below we will revisit this limit, by assuming the absence of anti-neutrinos in the decay of DM. For limits on the ν +ν flux from MeV-mass DM annihilation see [44]. resentative ideas), and therefore it is easy to anticipate that they can be used as tools for exploring DR.
In this paper we address several questions related to DR, WIMP direct detection experiments, and neutrino physics. Our primary goal is to assess-in broad classes of models-the sensitivity reach of direct detection experiments to DR. We shall concentrate on the eV-GeV energy range, and will mostly assume that a certain fraction of DM is prone to decay to DR. In the next section, we will outline several classes of possible interactions. In Sec. 3 we provide the computation of DR fluxes resulting from decaying DM, both in our halo and globally in the Universe. Section 4 considers the main phenomenological features of DR in several classes of generic models: DR in the form of the new population of neutrinos or neutrino-like light fermions, axion-like particles and light vectors ("dark photons"). We provide further discussion and conclusions in Section 5.

DR-SM INTERACTIONS
Among the possible sources of DR there can be processes involving the collision and decay of SM particles, and, under certain conditions, collisions, annihilation, and/or decay of massive DM particles. It is apparent that SM processes cannot create large amounts of DR, and in particular cannot saturate the last inequality in Eq. (2). The reason is that the only steady source of DR emitted globally in the Universe can come either from cosmic rays collisions with interstellar medium, or from the production of DR inside stars, including SN explosions. The energetics of these processes is subdominant to the energy locked inside ρ DM by many orders of magnitude.
On the contrary, very limited amount of information about the properties of DM provides some grounds for speculation that DM can indeed be a powerful source of DR. In the rest of this paper, we will concentrate on the decays of DM progenitor particles, giving rise to DR. We will call such a DM progenitor particle as X, and let κ be an (energy) fraction of DM that is allowed to decay. There are two generic types of DR we may consider. First, the DR can be the SM neutrinos: DM decay provides a late time source of energetic ν injection, and the SM-DR interaction in this case are the familiar neutrino weak interactions. Second, the quantum of DR could be a SM singlet, which we will label as χ; the pattern of the SM-χ interaction is more uncertain and diverse in this case. In the following we will discuss the general logical possibilities for these two cases.

A. SM neutrinos as dark radiation
In this scenario, we consider the possibility that DM gives rise to DR in the form of SM neutrinos. This can occur directly by X → ν(ν) + Y decays (where Y stands for the rest of the decay products). This can also occur in two steps: first through the decay of X to a nearly massless fermion χ (e.g., a sterile neutrino), that then oscillates into the SM neutrinos under certain conditions. Models with direct decay to neutrinos are free from potential constraints from N eff measurements, the ν-SM interactions are known, and the model is more minimal/simple. Then, on top of the constraints from dark matter experiments that limit the neutral current interactions of the DR neutrinos with nuclei, there will be additional constraints provided through weak interactions at neutrino detectors such as Super-Kamiokande (SK). The scenario with decays to neutrinos is also interesting in that it can be motivated by certain neutrino mass generation mechanisms, and connects to other aspects of neutrino physics and observables. Possibilities include: a vector or scalar boson DM particle X decaying to a pair of ν (ν), or fermionic DM X decaying to ν + Y where Y is a vector/scalar boson. Here we focus on the representative case where scalar DM X = φ decays to νν and/or νν.

Goldstone boson Majoron DM decaying to neutrinos
It has been considered that neutrino masses may arise from a global symmetry breaking (see e.g. [52,53] as well as [54] and references therein). The effective coupling of the Goldstone boson (Majoron) φ to neutrinos, resulting from the breaking of this global symmetry, is given by where L i stands for the SM lepton doublets, H is the SM Higgs boson, and Λ ij roughly corresponds to the global symmetry breaking scale. We have inserted the Higgs vacuum expectation value (vev) at the second step in the above equation. Note that here we retained the most generic flavor pattern in the coupling g ij which can account for neutrino mass and mixing. For simplification, we suppress the flavor indices of neutrinos and assume flavor-diagonal couplings in later discussion. In this case φ is real, and it decays to ν,ν symmetrically. The UV completion of the effective interaction in equation (3) can be the original Majoron model [52] where the authors assume a singlet Higgs field Φ responsible for the lepton number (L) symmetry breaking, and a Dirac pair of singlet fermions S L , S R which give rise to the familiar see-saw right-handed neutrino N once L-symmetry is broken. The Lagrangian relevant to m ν -generation is (we drop the lepton flavor indices) Upon symmetry breakings, this leads to two nonvanishing elements in the mass matrix of the system: m = y 1 H , M = y 2 Φ . Integrating out the heavy singlet fermions after symmetry breaking, one finds the prediction for neutrino masses and the Majoron-ν couplings as (assuming m M ): We can readily see how the parameters in the UV complete theory map to the effective coupling g in Eq.
(3): g y2m 2 ν y 2 1 H 2 . The mass of a Majoron as a pseudo-Goldstone has large uncertainty, generated by additional amounts of explicit global symmetry breaking (e.g. via a possible breaking of global symmetry by quantum gravity). There have been studies on the possibility of Majoron DM, and the focus has been on thermally produced Majoron, that typically needs to have mass O(keV) or less in order not to over-close the Universe [43,[55][56][57][58][59]. A Majoron with mass O(10) MeV or higher (which is our interest) can be relieved from this cosmological bound if it is produced non-thermally, through late decay of a massive particle (e.g. modulus) or late during inflation/reheating. One can check that realistic neutrino masses and a Majoron with lifetime around the age of the Universe can be simultaneously accommodated in this model, with perturbative couplings and Majoron masses in the range of our interest. In particular, the decay rate of φ in terms of effective parameter g or UV parameters can be written as follows: where t 0 = 4.3 · 10 17 s, the age of the universe, has been factored out for convenience.
A dedicated study of such a decaying Majoron DM and its relation to neutrino physics is an interesting topic on its own, but falls outside the scope of this paper.

Complex Scalar boson DM decaying to neutrinos
This scenario is inspired by the above Majoron model, yet instead we consider a complex scalar that carries lepton number but does not condense. Now the relevant interaction is: The UV completion of this model is similar to the Majoron model. But instead of Φ in the earlier model, we introduce a complex inert singlet, φ , which carries Lnumber but does not condense like Φ. Notice that the light scalar φ is not a Goldstone boson, so that its mass is not protected against radiative corrections. Thus, a sub-electroweak scale mass of φ may be associated with its own naturalness problem, which we do not address here. The relevant Lagrangian in this case has a similar form to that of the earlier model involving Φ: After the L -symmetry breaking, e.g. by Φ condensation, and after integrating out heavy states, we find: This model allows the possibility that there is an asymmetric abundance of φ vs. φ * , so that neutrinos instead of anti-neutrinos are produced as DR. Consequently, DR composed of ν rather thanν may not be subject to strong constraints imposed by SK on theν e flux just above the endpoint of the solar neutrino spectrum. The asymmetric DM type of φ requires an initial asymmetry generation (e.g. through the decay of a heavy Majorona RH neutrino) and efficient symmetric annihilation that depletes φ * . It is possible to realize such a model by extending the ingredients to include additional light (or massless) states enabling φ φ * annihilation.

DM decay giving rise to SM neutrinos through mass mixing
Another possibility of DM X decaying to SM neutrinos is through the mixing between SM neutrinos and light singlet sterile neutrinos that directly couple to DM. This falls into a similar category, as we will discuss in Sec. 2 B: X → χ + χ combined with the linear operator of mass mixing m χν χν SM . Depending on the χ − ν oscillation parameters, χ may or may not contribute to N eff . For example if the effective oscillation length between χ and ν is astronomically large, which can be achieved for nearly mass degenerate states χ and ν (and therefore very light χ), N eff does not provide an immediate constraint, as the ν → χ oscillation rate in the early Universe can be arbitrarily small. If in such a model the decays produce χ's of certain helicity, then the oscillations may result in the predominance of ν overν.

B. SM singlets as dark radiation
One can have several logical possibilities in this case: In general, such processes may occur with or without SM particles, or other members of dark sector particles Y in the final state, and the multiplicity of χ may vary. Also, the DR particles can be fermions or bosons. In other words, at this point there appears to be a great freedom in choosing models for X and χ, apart from γand cosmic-ray constraints if the decay is accompanied by energetic SM particles.
1. Singlet χ as a fermion: scattering signal Let us assume, for a moment, that DR χ is a fermion. Then, quite generically, the most important interactions of χ with the SM occur either at the linear or bi-linear order in χ, where O SM f (b) are generic fermionic or bosonic composite operators built from the SM fields; (χΓχ) stands for a generic bilinear operator in χ, and may contain a variety of currents, such asχγ µ χ,χχ etc. The overall transformation properties of O SM f (b) must be chosen such as to make these terms in the Lagrangian Lorentz invariant. To narrow down these many possibilities, we will consider where ν SM is some linear combination of SM neutrinos, and m χν is the effective mass parameter that will mix SM neutrinos and χ fermions. This mass parameter can be thought of as the low-energy limit of an operator respecting all symmetries of the SM. This way, χ will "participate" in the neutrino mass and mixing matrices. Among the bilinear interactions we will consider a subset of the most relevant ones for direct detection phenomenology, with two currents J EM and J B representing the SM electromagnetic current (defined without a factor of e), and the baryonic current of quarks. At the level of electrons, protons and neutrons, these currents can be approximated as where corrections related to the anomalous magnetic moments have been neglected. There are well-known UV completions of these operators where the exchange is mediated either by the dark photon (a vector particle V µ with kinetic mixing parameter to the photon field and Q χ g V µ (χγ µ χ) coupling), and/or by the baryonic vector (a new vector particle V B µ coupled to a baryon current via g V B ν J ν B ). In that case the effective couplings G V and G B can be expressed in terms of the more fundamental parameters as where Q χ is the charge of χ under the dark U (1) . While the dark photon models are fully UV consistent without further modifications, the baryonic model needs to be augmented at the weak scale to cancel the gauge anomaly in that sector. In the last decade, the phenomenological aspects of these models have been considered in some detail [9,[60][61][62][63][64][65][66]. Notice that recent works [67,68] have significantly advanced constraints on gauged baryon models, when V B µ is assumed to be light, and the gauge anomalies are canceled by a new set of fermions above the weak scale. Nevertheless such constraints do not directly apply, as (16) includes an extra parameter, Q χ , that can be large without violating perturbative unitarity (that requires Q χ g to be less than 4π). In this paper, however, we will stay on phenomenological grounds and treat G B as a free parameter. Also note that if mass of the mediator is zero in the G V coupling, χ will appear as fractionally charged fermion with effective EM charge of Q χ g . Milli-charged particles is an interesting and viable case of DR, as we will briefly discuss later.
While the choices given in (13) and (16) are far from being exhaustive, they are better motivated than many other ad-hoc models. Moreover, they are sufficient to capture the main phenomenological possibilities: ν − χ oscillation transitions and scattering of DR states on electrons and nuclei.

Singlet χ as a boson: absorption signal
If DR is represented by light bosonic particles, the possible choices of interactions are again plentiful. Natural possibilities of light bosons include dark photons and axion-like particles, both having masses protected by symmetries. These possibilities are also well-explored in the context of new directions for DM searches (see e.g. [69]). We first consider the decay of X to two dark photons, X → V V . This is a very economical model, where the interaction of dark photons with the SM is governed by the mixing parameter and the mass of the dark photon m V [70]. In vacuum, the dark photon couples to the electromagnetic current via V µ × (e J ν EM ). In medium, there is a well-known suppression of this interaction in the regime where m V is much smaller than the plasma frequency [71][72][73]. The other possibility of bosonic χ is an axion-like particle (ALP), a. The axion-like a can be produced through a scalar X DM decay X → aa, then a may interact with SM states via aF µνFµν (with the SM photon) or aēγ 5 e (with the electron). The main difference from fermionic DR is that bosons can be completely absorbed and converted to energy carried by the SM. The absorption of dark photons or ALP-type DM has been considered before in [45,49,74], and many features of these analyses can be generalized to ALP/dark photon DR detection.

C. Astrophysical and cosmological constraints on DR
In the cases where the DR considered here has a small or vanishing mass, there are many constraints coming from astrophysical and cosmological observations that one would need to take into account. Apart from the already mentioned CMB and BBN constraints, mainly in form of N eff , there are strong bounds imposed by stellar energy loss [75], which is particularly constraining for the case of X → aa, X → V V . In the case where the DR actually consists of massive (above O(MeV)) yet energetic particles such as boosted DM, BBN and CMB bounds do not automatically apply, although there can be model-dependent constraints, depending on how such massive DR interacts with the SM.
In several classes of models, DR that can induce observable effects for the direct detection and neutrino experiments without violating any of these indirect bounds.
Here we list some of the models that would clearly pass the indirect constraints: 1. Direct decay to neutrinos, X → νν for example, does not create any additional photons. Electronpositron pairs can be produced via X → ννe + e − but the rate will be suppressed by the weak interaction scale. Specifically, one may estimate that where the first factor reflects the additional phase space suppression, and the maximum available energy scale (m X ) is used to render the weak branching fraction dimensionless. Even for m X = 1 GeV, this is a ∼ 10 −13 level of suppression, which will not lead to additional constraints on the model from the global production of charged lepton cosmic rays.
2. X decay to χχ with bilinear interactions with the SM, Eq. (16). The same logic applies, and the probability of production of extra SM particles in a single decay process is suppressed by the phase space and in addition by G 2 V (B) m 4 X . (We will consider G V (B) in the ballpark of G F ). Here we assume, in addition, that the scale of the mediator mass is larger than the mass of the decaying particle, m V > m X , as otherwise a much more enhanced on-shell production of mediators may take place. With this extra condition, the indirect astrophysical constraints can also be completely avoided for sub-GeV X. These models are, however, constrained by N eff . The weak-strength (G V ∼ G F ) dark photon mediated interactions are capable to keep χ in a close thermal contact with electrons down to O(MeV) temperatures, creating a large positive contribution to N eff from the thermal bath of χ. Therefore, these models are disfavored unless G V G F , and then consequences of DR on direct detection are small. Models with χ coupled to the baryonic current fare better, as hadron number densities diminish rather quickly after the QCD phases transition. In this case, the meson-χ interaction could decouple relatively early, and G B on the order of G F is not manifestly excluded. For that reason we will concentrate on models with G B interaction, noting that the thermalization of new degrees of freedom via the baryon current interaction is a problem that deserves special separate study.
In some cases the indirect bounds from N eff can be circumvented at the expense of additional complication of models. For example, if the mass of χ is in the few MeV range, and, in addition, there is a light mediator particle in the same mass range, V , the self-annihilation χχ → V V with subsequent decay of V 's to the SM states may ensure that by the time of SM neutrino decoupling at MeV temperatures, all χ and V particles from the dark sector have annihilated and decayed, while X DM survives.
3. X decay to χχ, and the transfer of χ to neutrinos via χ → ν oscillations (a possibility that may be called "neutrino oscillation portal") can be made completely safe from cosmological and astrophysical bounds as m χν can be almost arbitrarily small. The smallness of this parameter does not imply a smallness of the χν mixing angle with the lowest neutrino mass eigenstate in the direct neutrino mass hierarchy picture, as there is no lower bound for the lightest m ν .
4. Finally, decays of X to nearly massless ALPs a, dark photons V , and light millicharged particles χ are not expected to be accompanied by a significant production of the SM states because of the extreme smallness of couplings between SM states and these types of dark sector particles.

MAXIMUM FLUXES OF DR FROM DECAYING DM
Any rate of detection of the relativistic background species χ (in this section χ stands for any DR particle, neutrino or singlet χ as discussed earlier) is controlled by the energy differential particle flux arriving at earth, dφ/dE χ . We consider two principal components of this flux, originating from DM decay within the galaxy ("gal") or from extragalactic distances ("eg"). 3 We will assume that only a fraction κ of the DM density decays (to be consistent with cosmological constraints given in [27] X is the lifetime of the progenitor X. Here, Ω dm = 0.1198h −2 is the CMB-inferred DM density parameter, and ρ c = 3M 2 P H 2 0 is the critical density today; H 0 = 100h km/s/Mpc with h = 0.6727 [1]. 4 Clearly, the fluxes attain their maximum when τ X is comparable to the age of the Universe, t 0 = 13.7 Gyr.
The decay of X injects particles χ with a spectrum, where N χ is the multiplicity of χ and Br χ is the branching ratio into χ in the decay with energy-differential width dΓ/dE χ . For a 2-body decay X → 2χ, the injection spectrum is a δ-function (N χ = 2), broadened only by the velocity dispersion of X prior to decay; for our purposes this is a negligible effect and in the following we take X at rest. The galactic particle flux at earth is found from the usual line of sight integral, Assuming no directional sensitivity, we average the Jfactor over all directions. Its value is relatively insensitive to the employed density profile and we use J 2.1 as obtained from a NFW profile; R sol = 8.33 kpc and ρ sol = 0.3 GeV/cm 3 are the distance to the galactic center and DM density at the position of the Earth.
For the extragalactic particle flux incident on Earth it is important to take the redshift of momentum into account. The flux, originating from the cosmological unclustered DM abundance, is given by the redshift integral, where the subscript "em" stands for the moment at emission. Equation (20)  diffusive extragalactic photon flux (in the limit of zero optic depth) when the mass of the daughter particle vanishes, m χ → 0. The only differences for m χ = 0 are contained in the properly blue-shifted energy and velocity at emission, obtained from the blue-shifted momentum, p em (z) = (1 + z) × p χ , where p 2 χ = E 2 χ − m 2 χ and p 2 em = E 2 em − m 2 χ . The cosmological redshift information in (20) is given through the Hubble expansion rate H(z) = H 0 (1 + z) 3 Ω m + Ω Λ and the cosmic time at redshift z, t(z). Since we are considering DM decays in the lowredshift Universe, it suffices to add contributions that originate from the matter-dominated era, z f < z eq ∼ 10 4 in (20). For a flat Universe one can then use the closed expression, with Ω m = 0.315 and Ω Λ = 1 − Ω m . At last, for a 2-body decay spectrum (18), the redshift integral in extragalactic flux can be evaluated directly to yield, where α ≡ p in /p χ ≥ 1. Note that the cosmic time t(z) is evaluated at redshift z = α − 1 in the exponential. Exemplary galactic and extragalactic fluxes are shown in Fig. 1 for τ X = 10 Gyr, m X = 50 MeV and κ = 0.1. A Gaussian of 5% width has been applied to the galactic flux for display. The total flux φ tot , integrated over the whole energy spectrum, varies over many orders of magnitude depending on the choice parameters. Nevertheless, one can estimate the maximum possible flux at κ ∼ 0.1, τ X = 10 Gyr, while taking m χ → 0, and keeping m X as a free parameter: Completely coincidentally, the value of the DR flux may become comparable to that of 8 B solar neutrinos at m X ∼ 10 MeV, and exceed diffuse SN neutrino flux by many orders of magnitude at m X ∼ 50 MeV. Fig. 2 demonstrates an example of total fluxes for varying lifetimes and X masses.

A. New population of SM neutrinos
If non-thermal DR radiation consists of SM neutrinos, we can predict their interaction rates in dark matter and neutrino detectors. The coherent nuclear recoil generated by the neutral current interaction is the easiest to treat, as it has no neutrino flavor or helicity dependence. The coherent neutrino-nucleus recoil cross section is given by, where Q W = Z(4s 2 w −1)+N is the weak charge of the nucleus. The scattering is essentially coherent in the number of neutrons N , owing to a cancellation in charge number Z since the weak angle is s 2 w 0.23. The degree of coherence is given by Helm form factor F (q) [79], evaluated at q = √ 2m N E R . Expression (24) is used by us to calculate the expected recoil signal in the DM direct detection experiments.
The coherent nuclear recoil is irrelevant for generating a signal in the most sensitive neutrino detectors. Instead, we must consider scattering on electrons (due to both, neutral and charged currents) and charged current scattering on nuclei. Moreover, there is a strong dependence of the expected signal on energy, flavor and helicity of DR neutrinos. In this paper we will assume dominance of neutrinos over antineutrinos in DR-a possibility discussed in Sec. 2 A 2-and, for simplicity, we will furthermore consider a flavor-universal content of DR at the interaction point.
From a few MeV to 1 GeV there are several approximate energy ranges that have significant differences with respect to the expected neutrino signal. These differences, to be discussed below, stem from a relatively large solar neutrino flux, and from different channels of neutrino interactions with electrons and nuclei.  [80]. Solar models predict [81] φ8 B = 5.46(4.50) × 10 6 cm −2 s −1 , depending on employed solar abundances [82] ( [83]). Taking the difference in the prediction as indicative of the size of the systematic error bar, one can conservatively, albeit somewhat ad-hoc, limit any additional flux of DR neutrinos in this energy range at nominal 90% C.L. as The main signature for neutrinos in this energy range is the scattering on electrons due to charged (CC) and neutral (NC) currents.
2. 16 MeV to 30 MeV. In this energy range the most important neutrino component isν e . The CC rate due to the inverse beta decay process on free protons, p +ν e → n + e + has a large cross section, and almost the entire anti-neutrino energy is communicated to the positron, E e + Eν − 1.8 MeV. Notice that there is no corresponding counterpart process for neutrinos, as n + ν e → p + e − is suppressed due to the binding energy of neutrons in carbon or oxygen nuclei. The resulting electron energy is at least ∼ 15 MeV lower than E ν , which puts these events in the energy range completely dominated by solar neutrinos. Thus, apart fromν e , in this energy range all neutrinos are detected due to their interaction with electrons inside the detector.
3. 30 MeV to 150 MeV. Starting at ∼30 MeV, reactions with neutrons inside nuclei (such as e.g. 16 O + ν e → 16 F + e) are no longer kinematically suppressed, and the energy of electrons in the final state (reduced relative to the incoming neutrino energy by ∼15 MeV or more) are above the energies of electrons created by solar neutrinos. In this energy range the CC cross sections of ν e andν e become similar and dominate over other forms of neutrino interactions. Above the muon threshold, the CC production of µ by ν µ also becomes possible, and it is constrained by the decay signature of stopped muons. In general, sub-Cerenkov energy muons irrespective of their origin represent a challenging background for searches of neutrinos of lower energy, as muon decays produce electrons that have sub-50 MeV energies. In our analysis, we account for both the CC interactions of ν e and the additional backgrounds introduced by muon decays.

Above 150
MeV. This is the energy range where the flux of atmospheric neutrinos is well measured, and found to agree with model predictions with ∼ 20% accuracy [38]. The main interactions forν e(µ) , ν e(µ) are the CC processes with nuclei. Therefore, using the results of [38], one can limit any extra DR neutrino component. Taking it conservatively, we demand that the DR flux shall not exceed ∼ 1/2 of the atmospheric neutrino flux, yielding dφ νe /dE < 0.08×E −2 GeV × cm −2 s −1 . (We stop our considerations of DR at about 1 GeV, and throughout this range a φ ν ∝ E −2 scaling holds reasonably well.) A constraint of similar strength applies to the muon neutrino flux above E ν = 300 MeV [38].
To summarize, the neutrino fluxes are least known directly in the energy regions 2 and 3. To treat constraints on DR with comparable energy range of neutrinos we require the expressions for the elastic and CC cross sections. The elastic scattering on electrons is given by where g L = −1 + 2s 2 w and g R = 2s 2 w . Indices e, µ, τ stand for the neutrino flavor dependence of the cross sections. For a flavor-universal composition, one should take σ el av = i=e,µ,τ σ el i /3. Because of the solar neutrino background, the most relevant quantity is the cross section that produces electrons above 15 MeV energy, The CC cross section ofν e can be calculated from first principles [84]. The same cannot be said about CC cross sections of ν e . Application to SK will require the cross sections for the 16 O + ν e → 16 F + e reaction as a function of energy. While it has not been directly measured over the whole energy range, one can use the results of theoretical calculations [85][86][87]. The cross section of a 30 MeV electron neutrino ν e on oxygen is 1.25×10 −42 cm 2 [86], rapidly decreasing below this energy. Various final states, such as different nuclear levels in 16 F and 12 B+α are possible, affecting the available measurable electron energy. We will treat this complication by assuming that the final state electron energy is, on average, shifted by ∼ 5 MeV below the threshold value, E e = E ν − 20 MeV and model the distribution by a Gaussian of 7% width. The modeling of this reaction is coarse-grained and can certainly be improved with a more dedicated study. Given the imperfect energy resolution of water Cerenkov detectors, we expect that the resulting limits will only depend mildly on our assumption.
The most important question to address now are the direct limits on DR neutrino fluxes that can be inferred from the SK data. As energy regions 1 and 4 are well understood in terms of neutrino fluxes, we need to determine the acceptable level of DR neutrino fluxes for regions 2 and 3. The experimental data relevant for this energy range are reported as a search for the diffuse supernova neutrino background [37]. While the SK collaboration applies their search to limitν e , owing to the fact that the CC cross section is largest, the same data can be used to limit other neutrino fluxes; relevant information can also be extracted from the high-end part of the 8 B solar neutrino spectrum [88]. In the 25-to-75 MeV range the background for the supernovaν e search is dominated by decay electrons that are produced from muon decays-themselves sourced from ν µ but undetected in their Cerenkov radiation-inside the SK volume, and by residual atmospheric electron neutrino CC events. A simultaneous fit to the shape of the signal plus background components allows the SK collaboration to extract a tight constraint on the flux ofν e .
We adopt a similar strategy and obtain a fit by taking into account the above sources of background with floating normalization (in addition to amplitude-fixed, smaller backgrounds that are inferred from "sidebands" to the SK search window; see Figs. 14-16 in [37]), together with the DR signal calculated via its elastic scattering on electrons and CC scattering of ν e on oxygen as outlined above; to the new component we apply the signal efficiency as reported in [37]. For the fit we minimize the likelihood ratio [89], where n i (µ i ) is the observed (expected) number of events in bin i and the sum runs over all bins of the SK search region (middle panels of Figs. 14-16 in [37]). In order to keep the scope of this investigation under control, in our numerical study we consider only the data from the SK-II period with 22.5 kton fiducial volume and 794 days of livetime; the inclusion of other SK runs is not expected to change the results qualitatively. An exemplary fit (pvalue 0.71), currently allowed by the data, is shown in Fig. 3. The red line is the DR signal, dominated by the CC reaction on oxygen.
The resulting likelihood fit of SK-backgrounds at each point in the (m X , τ x )-plane produces constraints on the admissible amount of DR. This constraint, at 95% C.L., is shown as the gray area labeled SK(dsnb) in Fig. 4, assuming that, close to its minimum, (28) follows a χ 2distribution. In addition, the gray areas SK(atm) and SK(sol) are excluded from atmospheric and solar flux measurements using (25) and reported fluxes in [88], respectively. The SK excluded regions supersede current constraints that are imposed by LUX [90] and XENON1T [91]. We derive the latter limits by computing the S1 scintillation signal from probability distribution functions that we derive with the statistical model described in [92]; see also the supplementary material of [93]. As input, the mean of the S1 signal is obtained from the light yield curve of Fig. 1 of [39] with a minimum nuclear recoil of 0.7 keV imposed. The employed overall light collection efficiencies are g 1 = 0.1 and 0.144, respectively. The generated signals are then constrained against data using the 'maximum gap method' [94] at 95% C.L.. To pick a specific example, we derive the current constraint on the neutrino flux originating from the decay of a DM particle of m X = 50 MeV and lifetime τ X t 0 : Notice that this constraint is more than two orders of magnitude more relaxed compared to the SK limit on a cosmicν e flux. Consequently, if this limit is saturated by DR, then the expected scattering rate inside the xenon-based direct detection experiments will be such as to mimic a DM particle with m DM ∼ 30 GeV and cross section of σ ∼ 10 −47 cm 2 , which is significantly above the traditionally derived "neutrino floor"; a more detailed numerical study on this point is in preparation. The SK constraints onν e can be significantly improved with the addition of Gd [30], as it would allow efficient detection of final state neutrons. This detector modification is unlikely to help strengthening constraints on other neutrino species, and therefore we project that SK constraints on ν DR are unlikely to be improved. On the other hand, the next generation of liquid xenon DM detectors is likely to reach the 10 −47 cm 2 level of sensitivity to DM, which would also make them efficient probes of neutrino DR. It has to be said that if the signal is detected at that level, one would have to check whether it is coming from scattering of nuclei on relativistic (DR) or non-relativistic (DM) species. This can e.g. be done, in principle, if another large and sensitive DM detector, based on a significantly lighter target nucleus, is built to complement xenon-based experiments. In particular, large scale argon-based detectors [95] may fill that niche. Constraints on SM neutrino DR, at 95% C.L.,in the assumption of X → νν decay, with equal weight for each neutrino flavor from SK measurements of atmospheric (atm) and solar (sol) fluxes as well as from searches for a diffuse SN neutrino background (dsnb). The limits supersede the current direct detection constraints derived from the final data set of LUX and from the initial data set of XENON1T; see main text for details.

B. New fermions interacting with SM via a dark force
Here we consider χ interacting through a vector portal that we take to be the baryonic current J µ B ( x), mediated by a massive vector V µ with mass m V . The incoming (semi-)relativistic particles χ will induce elastic scattering on nuclei in direct detection experiments. For that reason, we will generalize (24) to be valid both in the relativistic (E χ ≥ m χ ) and non-relativistic limits (E χ ≤ m χ ) alike. In order to compute the elastic recoil cross section dσ/dE R we note that the manifestly spin-independent (SI) part in the χ-nucleon matrix element p n |J µ B (0)| p n =ū p n Γ µ u pn is given by the vertex factor where F b (0) = 1 is the baryon number of the nucleon and in the following we can take it to unity; furthermore, we drop the term proportional to F m since it is suppressed. The direct detection recoil cross section on nuclei is then found to be The nuclear form factor F (q) is the same as above in (24). For m χ → 0, | p χ | → E χ and Eq. (31) reduces to the expression obtained in [96] using G B ≡ Q χ g 2 B /m 2 V and when the interaction obeys the contact limit q 2 m 2 V . The recoil rate in a direct detection experiment is then given by, where N T is the target number density per detector mass (for composite materials the rate is to be summed over the various target elements). The minimum χ energy E χ,min to produce a recoil E R is obtained from the corresponding minimum momentum | p χ,min | = E R m N /2. The predicted recoil rate in liquid xenon detectors is shown in Fig. 5. Current constraints from LUX [90] and XENON1T [91] on the model for G B = 10 G F is shown in Fig. 6. The solid (dotted) lines are for m χ = 0 (10 MeV); for details on the procedure that goes into the derivation of these limits cf. the preceding section.
To include constraints on baryonic-current coupled DR in χ from neutrino experiments, we have to consider the scattering of χ on nuclei that may lead to their recoil, break-up or excitation. We refer to previous studies [96], where it was shown that in the limit of small momentum transfer (in units of the nuclear size R N ), the interaction rate for inelastic processes is suppressed by (qR N ) 4 , which is a very small factor for q 100 MeV. This allows to tolerate large values of G B , without running into strong neutrino constraints. For this paper, we include constraints imposed on the model by scintillator-based neutrino detectors such as Borexino [97], that result from the elastic scattering on protons. The scattering of χ on protons, depending on E χ , can give a significant proton recoil that leads to energy deposition inside a liquid scintillator, and will give significant constraints on the higher end of the energy range of DR considered in this paper.
To calculate the proton recoil, we use (31) but take the nuclear form factor to be of the dipole form, F b (q 2 ) = (1+ q 2 /(0.71 GeV 2 )) −2 [98]. The treatment of the quenching factor in the recoil of protons is obtained following [99] and references therein; see also [96]. The gray region in Fig. 6 shows the resulting excluded region from measurements by Borexino, and it is derived from the fiducialized data shown in Fig. 3 of [100].
C. Discussion of ALPs and dark photons as DR Dark radiation in the form of ALPs or dark photons can appear as a result of X → aa or X → V V decays. The main question to investigate here is whether DM decay can provide a flux of ALPs or dark photons in the range where they can be detected while respecting limits from other, primarily astrophysical, constraints. In the case of ALPs, the most explored range is the keV frequency range, which for our scenario would imply m X ≥ O(keV), and the interaction with photons, L int = g aγγ a E · B. This energy range could make both, the dark matter experiments as well as axion helioscopes, sensitive to ALP DR. Previous studies have concluded that ALP DM with m a in the keV range is already severely restricted by the absence of a serious excess in X-rays (for a recent summary of X-ray constraints on decaying dark matter see e.g. [101]). The indirect constraints on ALP DM are typically much stronger than those provided by direct searches [45]. To remove constraints resulting from the decay of ALPs, one needs to require an additional hierarchy, m a m X . The easiest way to assess the detectability of DR in form of ALPs is to compare their maximal flux with the flux of solar axions. The spectral flux of the latter is given in [102], and integrating it over energy, one obtains φ solar a 10 12 (g aγγ × 10 10 GeV) 2 cm −2 s −1 . Comparing φ solar a to the maximum flux of ALPs attainable through DM decay, φ max tot in (23), we arrive at the maximum value of coupling g aγγ when the DR flux of ALPs have a chance of becoming larger than the solar flux, Current solar helioscopes utilize the a → γ conversion in the magnetic field to search for solar ALPs. The same techniques can be used to search for ALPs forming DR. (For some values of parameters, the galactic component of DR dominates, and one should expect an enhancement of the conversion in the direction to the galactic center.) It is easy to see that in the keV range of frequencies, the benchmark value g aγγ = 10 −11 GeV −1 in Eq. (33) is outside the reach of the current generation keV-range ALP detector CAST [102], but may be amenable to searches with the next generation ALP telescopes, such as IAXO [103]. A similar conclusion can be reached for ALPs coupled to the electron spin via g aee a ×ēγ 5 e. The expected solar flux of ALPs is at the level of ∼ (g aee × 10 13 ) 2 × 10 9 cm −2 s −1 [104], which is again somewhat larger than the maximum attainable flux for DR with a keV scale progenitor X, unless g aee is below 10 −13 . The current sensitivity of dark matter experiments to solar axions is at the level of g aee ∼ 8 × 10 −12 [105], and, therefore, only significant improvements in the sensitivity of large-scale dark matter experiments could render a hypothetical ALP component of DR detectable.
We now turn to the case of dark radiation in the form of dark photons. The main difference with the ALP case is that for small mass of dark photons, the solar flux decouples as φ sol V ∝ m 2 V [72]. (See Ref. [106] for a detailed calculation of the solar energy loss to dark photons.) On the other hand, the production of DR dark photons may not need to be suppressed by small mass m V in the same limit, and therefore the flux of DR can be parametrically larger than the solar flux of dark photons. The analysis of the absorption of dark photon dark radiation is very similar to the analysis performed for dark photon DM [49,74]. The conclusion of these studies is that in some corners of mass-mixing angle parameter space, {m V , }, the direct detection experiments have sensitivity to dark photons beyond the astrophysical constraints. Unfortunately, our analysis shows that the dark photon dark radiation is currently not constrained by direct detection experiments. The qualitative reason for that is as follows: when the medium effects can be neglected (m V > 1 eV) the absorption rate per atom scales as 2 ×n V ×σ photo ×c, where n V is the number density of dark photons (in the form of DM or DR), σ photo is the photo-ionization cross section, and c is the speed of light. The rate is approximately independent on the velocity of dark photons when their total energy is fixed. For dark photon DM the number density n V is typically larger than for dark photon DR, resulting in a smaller ionization rate for the latter. Taking a representative point on the parameter space, = 10 −11 , m V = 1 eV, and κ = 0.1, m X = 50 eV, and τ X = t 0 , we find that the ionization rate for the XENON10 experiment is not exceeding a few times 10 −7 /kg/day, which is much smaller than the current sensitivity.
Finally, if the dark photon mass is very small or zero, it can mediate the interaction between charged particles of the SM and the dark sector χ particles when they are charged under the dark photon gauge symmetry. Such objects, generically called milli-charged particles, have been extensively studied in the literature with the latest constraints compiled in [107]. While the range of small masses is again generically very constrained by the combination of cosmology and astrophysics, the GeV and heavier range χ are allowed to have charge up to (10 −2 − 10 −1 ) × e. If the DM particle X decays to a pair such millicharged particles, one should expect a variety of new effects associated with scattering of χ inside dark matter and neutrino detectors.

CONCLUSIONS
We have considered a hypothetical possibility that along with non-relativistic DM, some (semi-)relativistic particles form a cosmic dark radiation (DR) background that may have a noticeable interaction rate with the SM particles. Such DR can be a non-thermal component of energetic neutrinos, or SM singlet particles. The most efficient mechanism for populating DR radiation states is the decay of DM, and if it happens at late redshift, a significant fraction of DM is allowed to decay to DR. Adopting this framework, we have derived the energy spectrum of DR that includes the galactic and a global cosmological component, as a function of the progenitor particle's lifetime, mass and abundance.
To narrow the discussion we have concentrated on the sub-GeV range for the energy of DR (or, equivalently, the mass of decaying DM particles). This range is the most relevant one for a potential signal in the detectors that are built to search for WIMP DM recoils. Therefore, they could also register energy and momentum transfer that is communicated to nuclei and electrons via their interaction with DR. For the masses of decaying particles in the tens of MeV range, the resulting fluxes of DR particles may reach 10 5 − 10 6 cm −2 s −1 , which is a fairly significant flux, exceeding atmospheric and expected diffuse SN neutrino fluxes by many orders of magnitude. Correspondingly, we find that both, neutrino and DM direct detection experiments are sensitive to the weakscale interaction of DR with SM fermions (nucleons and electrons).
For DR in the form of the SM neutrinos, we find that the Super-Kamiokande (SK) experiment provides the dominant constraints. The strongest constraints stem from limits onν e fluxes. If, however, DM decays preferentially to neutrinos, rather than antineutrinos, the constraints become much milder, and current SK data can tolerate much larger fluxes of DR neutrinos. If one saturates our derived limits on DR neutrino fluxes, DM direct detection experiments need to improve their sensitivity by approximately two orders of magnitude to become competitive with SK. Another way of phrasing the same finding is to say that the signal from DR-induced neutrino-nucleus scattering could be significantly above the normally expected "neutrino floor" for DM direct detection experiments. In particular, we find that for the energy of DR particles in the ∼ 25 MeV range, the gain over the nominal neutrino floor can be very significant, providing a potential signal that would compare to a one from a ∼ 30 GeV WIMP-type dark matter particle with a scattering cross section of 10 −47 cm 2 .
Among the most interesting cases for DR in the New Physics sector are the dark radiation models that interact with nuclei via a baryonic current, which is a possibility that is least constrained by neutrino experiments. Here, the existing direct detection experiments provide dominant bounds for the same DR energy range. Dedicated studies with upcoming neutrino experiments such as SNO+, JUNO, and DUNE/LBNF can potentially be complementary probes for these models.
To conclude, the co-existence of DM and interacting DR is a generic possibility that draws an analogy with the structure of the SM, where both massive particles (atoms) and radiation (photons, neutrinos) are present. This work demonstrates that this broad new class of physics can be probed with experiments originally designed to search for dark matter and study neutrino interactions.