Neutrino Discovery Limit of Dark Matter Direct Detection Experiments in the Presence of Non-Standard Interactions

The detection of coherent neutrino-nucleus scattering by the COHERENT collaboration has set on quantitative grounds the existence of an irreducible neutrino background in direct detection searches of Weakly Interacting Massive Dark Matter candidates. This background leads to an ultimate discovery limit for these experiments: a minimum Dark Matter interaction cross section below which events produced by the coherent neutrino scattering will mimic the Dark Matter signal, the so-called \emph{neutrino floor}. In this work we study the modification of such neutrino floor induced by non-standard neutrino interactions within their presently allowed values by the global analysis of oscillation and COHERENT data. By using the full likelihood information of such global analysis we consistently account for the correlated effects of non-standard neutrino interactions both in the neutrino propagation in matter and in its interaction in the detector. We quantify their impact on the neutrino floor for five future experiments: DARWIN (Xe), ARGO (Ar), Super-CDMS HV (Ge and Si) and CRESST phase III (CaWO$_4$). Quantitatively, we find that non-standard neutrino interactions allowed at the $3\sigma$ level can result in an increase of the neutrino floor of up to a factor $\sim 5$ with respect to the Standard Model expectations and impact the expected sensitivities of the ARGO, CRESST phase III and DARWIN experiments.

Abstract: The detection of coherent neutrino-nucleus scattering by the COHERENT collaboration has set on quantitative grounds the existence of an irreducible neutrino background in direct detection searches of Weakly Interacting Massive Dark Matter candidates. This background leads to an ultimate discovery limit for these experiments: a minimum Dark Matter interaction cross section below which events produced by the coherent neutrino scattering will mimic the Dark Matter signal, the so-called neutrino floor. In this work we study the modification of such neutrino floor induced by non-standard neutrino interactions within their presently allowed values by the global analysis of oscillation and COHERENT data. By using the full likelihood information of such global analysis we consistently account for the correlated effects of non-standard neutrino interactions both in the neutrino propagation in matter and in its interaction in the detector. We quantify their impact on the neutrino floor for five future experiments: DARWIN (Xe), ARGO (Ar), Super-CDMS HV (Ge and Si) and CRESST phase III (CaWO 4 ). Quantitatively, we find that non-standard neutrino interactions allowed at the 3σ level can result in an increase of the neutrino floor of up to a factor ∼ 5 with respect to the Standard Model expectations and impact the expected sensitivities of the ARGO, CRESST phase III and DARWIN experiments.

Introduction
The existence of Dark Matter (DM) in the Universe is considered one of the pillars of physics Beyond the Standard Model (BSM). However, almost a century since it was suggested we remain ignorant of what it is made of. Nature has not been kind enough so far to provide us with an unambiguous way to probe its character by means of its interaction with a human made detector. This is so despite the continuous efforts of scientific and engineering ingenuity which has resulted into experiments which are sensitive to weaker and weaker DM interactions and to a larger range of the DM candidate masses, in particular, for the well motivated DM models consisting of Weakly Interacting Massive Particles (WIMPs) [1].
Neutrinos, on the other hand, have given us the first direct evidence of BSM physics in the form of flavor transitions which were first observed in natural neutrino fluxes produced in the Sun and in the Earth atmosphere and have now been confirmed and precisely measured in a variety of oscillation experiments [2]. As usual in particle physics the signals of today are the backgrounds of tomorrow, and the possibility of the coherent scattering of these natural neutrinos with nucleus (CνNS) [3] of the target material giving a nuclear recoil, similar to what is expected of a DM particle, will eventually constitute an irreducible background for the WIMP discovery [4] in the next generation of DM direct detection experiments [5][6][7][8][9][10][11]. These experiments are expected to have sensitivity to significantly lower WIMP-nucleon cross sections thanks to a lower recoil energy threshold and/or a larger target mass. This, in turn, will place these detectors closer to reach such discovery limit also known as the neutrino floor. Thus, the first observation of CνNS by the COHERENT collaboration [12] with rates consistent with the SM expectation bears important implications for them.
In this paper we want to quantify what is the range of variation of the neutrino floor for such next generation experiments which is still allowed by the present experimental constraints on the relevant neutrino interactions. One can phenomenologically study generic contributions that affect CνNS, and therefore can modify the neutrino floor, by considering a model-independent approach where flavor dependent four-fermion effective operators of neutrinos and quarks are allowed. These operators will induce neutral current (NC) nonstandard neutrino interactions (NSI) with the target nucleus at the detector, and also modify in a non-trivial manner the matter potential relevant to the flavor evolution of the solar and atmospheric neutrinos. For this reason not only COHERENT scattering data but also results from oscillation experiments are relevant in the determination of their allowed values.
There are several studies in the literature of how the expected rate of neutrino background in direct DM detection experiments could be altered by non-standard flavor-independent contributions [13,14] as well as flavor-dependent NSI's [15,16]. Here we extend those works with a detailed and statistically consistent evaluation of the impact of NSI's on the neutrino floor at a variety of next generation DM detection experiments. To do so we make use of the full likelihood information of the global analysis of oscillation and COHER-ENT data to consistently account for the correlated effects of the NSI's both for neutrino propagation in matter and in its interaction in the detector. In particular, we quantify their impact on the neutrino floor for five future experiments: four with monoatomic detector targets (Xe used by DARWIN, Ar by ARGO, Ge and Si by Super-CDMS HV) and one operating with a molecular target (CaWO 4 for CRESST phase III).
This paper is organized as follows. In Sec. 2 we introduce the relevant NC-NSI Lagrangian and summarize the results on the constraints on the NC-NSI coefficients from the global analysis of neutrino oscillation and COHERENT scattering data. In Sec. 3, we use the full likelihood of such global analysis to quantify the allowed modification to the differential recoil rate of neutrino events at direct detection experiments by the presence of these NC-NSI. After describing briefly the full statistical analysis used to obtain the DM discovery limit, in Sec. 4 we present our results in terms of the range of variation of the neutrino floor induced by the NC-NSI interactions at the aforesaid future detectors using different target materials. Finally, in Sec. 5 we draw our final conclusions.

NC-NSI Lagrangian and constraints
In this work we consider NC-NSI affecting neutral-current processes relevant to neutrino coherent scattering in DM experiments. Generically NC-NSI interactions can be parametrized in the form: where G F is the Fermi constant, α, β are flavor indices, P ≡ P L , P R and f is a Standard Model (SM) charged fermion. In this notation, f,P αβ parametrizes the strength of the new interaction with respect to the Fermi constant, f,P αβ ∼ O(G X /G F ). First we notice that if all possible operators in Eq. (2.1) are added to the SM Lagrangian, the Hamiltonian of the system which governs neutrino oscillations in presence of matter is modified as 2) where U vac is the 3-lepton mixing matrix in vacuum [17][18][19] and N e (x) is the electron number density as a function of the distance traveled by the neutrino in matter. For antineutrinos Hν = (H vac −H mat ) * . In Eq. (2.2) the generalized matter potential depends on the effective NC-NSI parameters, αβ , defined as Note that the sum only extends to those fermions present in the background medium (up quarks, down quarks and electrons), and is the average ratio of the number density of the fermion f to the number density of electrons along the neutrino propagation path. In the Earth, the ratios Y f are constant to very good approximation, while for solar neutrinos they depend on the distance to the center of the Sun. The presence of NC-NSI with electrons, f = e, would affect not only neutrino propagation in matter, as described in Eq. (2.2), but also the neutrino-electron cross section in experiments such as SK, Borexino and reactor experiments; hence, for the sake of simplicity, in what follows we consider only NC-NSI with either up quarks (f = u) or for down quarks (f = d). Also, it should be noted that only the vector-like NC-NSI combination ( V = L + R ) contribute to the matter potential in neutrino oscillations. Those are precisely the relevant couplings for CνNS in the approximation E R E ν , where E R is the nuclear recoil energy and E ν the incoming neutrino energy.
In general, the matter potential in Eq. (2.2) contains a total of 9 additional parameters per f : 3 diagonal real parameters and 3 off-diagonal complex parameters (i.e., 3 additional moduli and 3 complex phases). However, the evolution of the system given by the Hamiltonian in Eq. (2.2) is invariant up to a constant, so that oscillation experiments are only sensitive to the differences between the diagonal terms in the matter potential.
As a consequence of the CPT symmetry, neutrino evolution is invariant if the Hamiltonian in Eq. (2.2) is transformed as H ν → −(H ν ) * , see [20,21] for a discussion in the context of NC-NSI. This transformation can be realized by changing the oscillation parameters as 4) and simultaneously transforming the NC-NSI parameters as see Refs. [21][22][23] for details. In Eq. (2.4) δ CP is the leptonic Dirac CP phase, and we are using here the parametrization conventions from Ref. [23]. This degeneracy allows for a change in the octant of θ 12 , which leads to the appearance of the so-called LMA-D solution in the solar neutrino analysis [24], when combined with large diagonal NC-NSI coefficients ee − µµ ∼ −2, and it also implies an ambiguity in the neutrino mass ordering. Being an intrinsic degeneracy of the full OSC+NSI Hamiltonian for neutrino propagation it cannot be resolved by neutrino oscillation experiments alone.
The presence of the NC-NSI's also affects the inelastic interaction of neutrinos and therefore neutrino scattering experiments are sensitive to some combinations of f,P αβ , depending on whether the scattering takes place with nuclei or electrons, the number of protons and neutrons in the target nuclei, and other factors. Thus altogether current experimental constraints on vector-like NC-NSI parameters include those obtained from a global fit to oscillation data [21] as well as those obtained from results from neutrino scattering data. For oscillation constraints on NC-NSI parameters we refer to the comprehensive global fit in the framework of 3ν oscillation plus NC-NSI with up and down quarks performed in [21] which we briefly summarize here for completeness. All oscillation experiments except SNO are only sensitive to vector NC-NSI via matter effects as described above. There is some sensitivity of SNO to axial couplings in their NC data, so for the sake of simplicity the analysis in Ref. [21] and all combinations that we will present in what follows are made under the assumption of purely vector-like NC-NSI. The fit includes data sets from KamLAND reactor experiment [25] and solar neutrino data from Chlorine [26], Gallex/GNO [27], SAGE [28], Super-Kamiokande [29][30][31][32][33] Borexino [34,35] and SNO [36][37][38][39], together with atmospheric neutrino results from Super-Kamiokande phases 1-4 [40], long-baseline results from MINOS [41,42] and T2K [43], and reactor results from CHOOZ [44], Palo Verde [45], Double CHOOZ [46], Daya Bay [47] and RENO [48], together with reactor short baseline flux determination from Bugey [49,50], ROVNO [51,52], Krasnoyarsk [53,54], ILL [55], Gösgen [56], and SRP [57].
In principle, the analysis depends on the six 3ν oscillations parameters plus eight NC-NSI parameters per f target. To keep the fit manageable in Ref. [21] only real NC-NSI were considered and ∆m 2 21 effects were neglected in the analysis of atmospheric and long-baseline experiments. This renders the analysis independent of the CP phase in the leptonic mixing matrix. Furthermore in Ref. [58] it was shown that strong cancellations in the oscillation of atmospheric neutrinos occur when two eigenvalues of H mat are equal, and it is for such case that the weakest constraints are placed.
For what concerns bounds from non-oscillation experiments, generically NC-NSI parameters in gauge invariant models of new physics at high energies are expected to be subject to tight constraints from charged lepton observables [59,60] and also directly from collider data (see for example [61]). For those constructions the bounds on the NC-NSI of neutrinos are too strong to lead to any observable effect in direct dark matter detection experiments. However, more recently it has been argued that viable gauge models with light mediators (i.e., below the electro-weak scale) may lead to observable effects in oscillations without entering in conflict with other bounds [62][63][64] (see also Ref. [65] for a discussion). Furthermore, for very light mediators, bounds from high-energy neutrino scattering experi-ments such as CHARM [66] and NuTeV [67] do not apply. Conversely, even for models with such very light mediators bounds on the vector-like NC-NSI can be imposed by the direct observation of coherent neutrino-nucleus scattering, i.e., the same reaction which would give rise the neutrino floor in the DM direct detection experiments. The modification to the CνNS cross section is given by [68] where F(E R ) is the Helm form factor [69], m N is the nuclear mass, and E ν and E R are the neutrino and recoil energies, respectively. The factor Q α NSI depends on the f,V αβ parameters In this respect, the most relevant results are those from the COHERENT collaboration which recently reported the observation of coherent neutrino-nucleus scattering at 6.7σ [12] with neutrinos produced from pion decay at rest on a 14.6 kg CsI[Na] detector target. The observed interaction rate is in good agreement with the SM prediction and therefore can be used to constrain NC-NSI. In Ref. [70] the COHERENT results were combined with the neutrino oscillation constraints to produce model-independent bounds on vector-like NC-NSI parameters which are relevant to the evaluation of the neutrino floor for DM detection. In particular, it was shown that after the inclusion of the COHERENT data the global analysis excludes the intrinsic degeneracy in the NSI+OSC Hamiltonian described above (i.e., the LMA-D solution) at 3.1σ (3.6σ) CL for NC-NSI with up (down) quarks. In addition, the combination of oscillation and COHERENT allows to derive competitive constraints on all NC-NSI parameters. For illustration we show in Fig. 1 the results from that combined fit by plotting αβ , for f = u (upper panels) and f = d (lower panels) after marginalization over the undisplayed oscillation and NC-NSI parameters in each panel. 1 As seen in the figure the combination of oscillation and COHERENT results is able to impose constraints on all the relevant vector-like NC-NSI coefficients even when they are varied simultaneously in the analysis. In this respect, it is important to notice that, in order to consistently obtain the allowed effects on the neutrino floor at a given CL as presented in the following section, we use the full multidimensional ∆χ 2 dependence on the NC-NSI and oscillation parameters so to correctly account for the correlations among their allowed ranges.  Figure 1. ∆χ 2 as a function of the NC-NSI parameters f,V αβ , for a global fit to oscillation experiments and COHERENT data (see text for details).

Neutrino rate at direct detection experiments
Let us start by evaluating the impact of NC-NSI on the differential recoil rate of solar and atmospheric neutrinos at direct detection searches. In doing so one has to bear in mind that, unlike for SM ν interactions [4] and for the simplified BSM models considered in Ref. [14], NC-NSI are flavor dependent and therefore in their presence one must take into consideration the flavor of the neutrino fluxes arriving at the detector. This is incorporated by introducing the neutrino transition probability between the source and the detector, so the differential recoil rate takes the form: where N T corresponds to the total number of targets in the detector, dΦ dE ν | να is the incoming flux of neutrinos with flavor ν α , P (ν α → ν β , E ν ) is the transition probability between the flavors ν α → ν β , dσ ν (ν β ) dE R | NSI is the flavor dependent CνNS cross section given in Eq. (2.6) and E ν min is the minimum neutrino energy that a neutrino must have to create a recoil with energy E R , that is For simplicity, we will only consider here the solar and atmospheric contributions to the neutrino flux. The Diffuse Supernova Neutrino background [72], expected but not yet observed, could affect the discovery limit for very large exposures [73,74]. Other neutrino fluxes, such as geoneutrinos [4] and reactor antineutrinos [75], depend on the location of the detector on the Earth, but in any case their intensity is very small at practically all the important sites, thus making their contribution negligible for the facilities and exposures analyzed here. The flavor transition probabilities for solar ν e and atmospheric ν e and ν µ (and their antineutrinos) in the presence of NC-NSI are obtained by solving the evolution equation: with H ν given in Eq. (2.2). Solar neutrinos are electron neutrinos produced in the Sun core with energy up to ∼ 20 MeV. They are usually labeled as pp, pep, hep, 7 Be, 8 B, 13 N, 15 O, and 17 F neutrinos depending on the reaction rate in which they are produced. In what follows we will consider the fluxes as predicted by the B16-GS98 Solar Model [76]. Its flavor transition probabilities can be obtained in very good approximation in the one mass dominance limit, ∆m 2 31 → ∞ (which effectively means that generically G F N e Y f f,V αβ ∆m 2 31 /E ν ). In this approximation the survival probability P ee can be written as [77,78] P ee = c 4 13 P eff + s 4 13 , (3.6) where for δ CP we have used the convention of Ref. [23]. In this scheme, the coefficients f Atmospheric neutrinos are electron and muon neutrinos and antineutrinos produced by the interaction of cosmic rays with the atmosphere and can extend to very large energies. In what follows we consider the fluxes obtained by the FLUKA collaboration [79]. Only the low energy (sub-GeV) part of their spectrum can be relevant as background in DM direct detection experiments [4]. Transition probabilities for these low energy neutrinos present some unique features [80,81] but generically Earth matter effect are expected to have limited impact in their evolution. Nevertheless, for the sake of consistency, we obtain the flavor transition probabilities for the atmospheric neutrinos by solving the 3ν evolution equation (Eq. (3.3)) using the Earth matter density distribution described by the PREM model [82]. Since direct detection experiments do not have directionality, one has to consider neutrino trajectories arriving in the detector from any direction. Both the length of the trajectory and the amount of Earth matter traversed will be different for neutrinos arriving from different zenith angles, which makes their flavor probability at the detector directiondependent. The transition probability included in the differential recoil rate, Eq. (3.1), is then the angular average In Fig. 2 we plot the differential recoil rate of solar and atmospheric neutrinos as a function of the recoil energy for a monoatomic target (Xe, upper panels) and for a molecular one (CaWO 4 , lower panels). We only show the effects of NC-NSI for f = u as the cases for f = d are very similar. For illustration, we label in the figure the dominant component of the neutrino flux contributing at the different targets for each recoil energy range. The shown ranges are obtained by scanning the multidimensional space of NC-NSI and oscillation parameters allowed by the global analysis of oscillation+COHERENT data at 1σ, 2σ, and 3σ (∆χ 2 OSC+COH ≤ 1, 4,9). For each point in the space of the allowed parameters we compute the predicted rate consistently including the NC-NSI coefficients both in the oscillation probabilities and in the interaction cross section.
We notice first that the prediction for best fit parameters is off the SM expectation. As mentioned in Sec. 2 this is driven by the 2σ tension between the determination of ∆m 2 21 from KamLAND and solar neutrino experiments. From the figure we also see that the allowed modification at 1σ, 2σ and 3σ can result both in an increase or a reduction of the expected rate. Quantitatively we find that with NC-NSI's allowed at 3σ by the oscilla-tion+COHERENT data the neutrino background rate can vary from a factor ∼ 6-7 larger to a factor ∼ 0.2 smaller than what is predicted by the SM interactions, the larger enhancement happening for a CaWO 4 detector. Moreover, as seen in the figure, the modification of the rate with respect to the SM prediction depends on which nuclear reaction dominates the solar neutrino flux at the given energy. This happens because of the NC-NSI in the matter potential as well as the different distributions of neutrino production points for the various nuclear reactions, which implies that the modification of the transition probabilities depends on the neutrino energy. However, altogether this is a relatively small effect. For solar neutrinos, the inclusion of the NC-NSI in the matter potential leads to correction of up to ∼ 1% (∼ 12%) in the rate of pp ( 8 B) neutrinos, while for atmospheric neutrinos the   Figure 2. Expected differential recoil event rate (left) and the fraction of the NSI contribution compared to the SM one (right) as a function of the recoil energy E R for the solar and atmospheric neutrino fluxes. In the upper (lower) panels the detector target is Xe (CaWO 4 ). The full line corresponds to the prediction for the best fit NC-NSI and oscillation parameters of the analysis of oscillation+COHERENT data, and the shaded regions bordered by dashed, dot-dashed, and dotted lines correspond to predictions with NSI and oscillation parameters with ∆χ 2 OSC+COH ≤ 1, 4, 9, respectively. difference induced by including NC-NSI in the propagation is of order ∼ 1%. Finally, note that the peak structure that comes up at E R ≈ 50 keV for CaWO 4 is a form factor feature. Atmospheric neutrinos are relevant for recoil energies above O(30 keV).

Results: NC-NSI effects on the direct detection discovery limit
A WIMP direct detection experiment is based on a simple principle: if a WIMP interacts with the quarks of a nucleus, it should be possible to detect such interaction by measuring the nuclear recoil produced on a target material. The negative results obtained so far have propelled several proposals for future experiments [5][6][7][8][9][10], which will have larger exposures and will be sensitive to smaller WIMP masses. It is for such future experiments that the solar and atmospheric neutrino background can represent a limitation to sensitivity.
Bounds on the WIMP mass-cross section parameter space accessible by a given experiment can be calculated using various statistical methods. In this work we follow the approach of Refs. [4, 73, 83? , 84] which makes use of the likelihood ratio test statistics to calculate the so-called discovery limit, which is defined as the minimum cross section for a given WIMP mass, m χ , for which a 3σ discovery is possible in 90% of the hypothetical experiments. In brief, we define a binned likelihood test statistics with n b = 50 bins in the interval [E th , 100 keV], where E th is the energy threshold of the experiment: with P being the Poisson distribution function for measuring a N i obs number of events in bin i from the expected number of WIMP+neutrino events N i χ + nν j=1 N i ν (φ j ν ). The sum over neutrino species j includes all solar and atmospheric components. The expected number of neutrino events in bin i from neutrino flux component φ j ν is trivially computed by integrating the corresponding differential recoil rate in the bin, while the expected number of WIMP events N i χ is 3) The WIMP differential recoil rate is [69] dR where ρ 0 = 0.3 GeV/c 2 /cm 3 is the local DM density and µ n = m n m χ /(m n + m χ ) is the reduced mass of the WIMP-nucleon system [69,85]. Furthermore, we assume a Maxwell-Boltzmann distribution function for the WIMP velocity distribution in the Earth's frame of reference f (v), with the escape velocity v esc = 544 km/s, and the Earth's velocity v lab = 232 km/s [69]. In Eqs. (4.2)-(4.3) ε(E R ) is the detector efficiency function. The terms L(φ j ν ) in Eq. (4.1) correspond to prior factors included to parametrize the theoretical uncertainties on the normalization of each neutrino flux component. We take them to be Gaussian distributed with uncertainties σ j given by the B16-GS98 Solar Model (6% for pp, 10% for pep, 30% hep, 6% for 7 Be, 12% for 8 B, 15% for 13 N, 17% for 15 O, and 20% for 17 F), and the atmospheric flux calculations which we conservatively take to be 20% [86]. These neutrino flux normalizations will be treated here as nuisance parameters, while the energy spectrum of each component is assumed to be known. 2 For the sake of simplicity we do not include in our likelihood astrophysical uncertainties affecting the WIMP signal [74].
The test between the neutrino-only hypothesis, H 0 , and the neutrino+WIMP hypothesis, H 1 , is performed by using the likelihood ratio for a fixed WIMP mass whereφ ν andσ χn are the neutrino fluxes and WIMP-nucleon cross section values which maximize the likelihood L(σ χn , m χ ;φ ν ) whileφ ν is obtained by maximizing the likelihood function in the case of the null hypothesis, L(σ χn = 0, m χ ;φ ν ). The significance of the WIMP signal (or, equivalently, of the disagreement between H 0 and H 1 ) is Z = −2 ln λ(0). We construct its probability distribution function, p(Z|H 0 ), by MC simulating 2000 experiments, creating for each one a spectra for both WIMP and neutrinos and computing their corresponding value of Z. To determine the discovery limit we need to find the smallest cross section for which 90% of the experiments have evidence larger than 3σ, this is, the minimum value ofσ χn for which Clearly for an experiment with a given target, exposure and threshold the presence of NC-NSI will affect the discovery limit as it modifies the background neutrino event rate. To illustrate the maximal possible size of the modification for a wide range of WIMP masses we start by considering an idealized experiment with an extremely low threshold (E th = 0.01 eV), a large exposure (10 3 ton-year) and a 100% efficiency (ε(E R ) = 1). We plot the results in Fig. 3 for the same two targets considered in Fig. 2, i.e., a monoatomic target of Xe and a composed target made of calcium tungstate (CaWO 4 ). For the sake of concreteness we consider NC-NSI with f = u quarks but similar results hold for NC-NSI with f = d.
As seen in Fig. 3 the discovery limit curves exhibit the well known peaks [4, 73, 84? ] corresponding to the region in the parameter space which is significantly mimicked by the neutrino flux component as labeled. For instance, the peak near a WIMP mass of 6 GeV is related to the background created by 8 B solar neutrinos. Again, the shown ranges are obtained by scanning the multidimensional space of NC-NSI and oscillation parameters allowed by the global analysis of oscillation+COHERENT data at 1σ, 2σ and 3σ (∆χ 2 OSC+COH ≤ 1, 4,9). For each point in the allowed parameter space we compute   the predicted neutrino background rates in each bin, consistently including the NC-NSI coefficients both in the oscillation probabilities and in the interaction cross section. With such rates we obtain the corresponding discovery limit. As expected the inclusion of NC-NSI does not change the position of the peaks since they do not modify the E ν dependence of the CνNS cross section. Also, as mentioned in the previous section the neutrino floor predicted at the best fit point does not exactly coincide with the SM prediction because of the preference of oscillation data for a small but non-zero NC-NSI coefficients, which helps reducing the tension between solar and KamLAND data in the determination of ∆m 2 21 . As seen in the figure NSI-NC can lead to an enhancement or a suppression of the discovery limit. Quantitatively we find that within the 3σ allowed region the neutrino floor can be modified between a factor ∼ 3 ( 13 N neutrinos) to ∼ 5 (hep neutrinos) with respect to the expectation from the SM. Finally, comparing the two panels in Fig. 3 we see that the main difference between the two targets appears in the region of WIMP masses around ∼ 10 GeV, and is a consequence of the different recoil energy behavior of the neutrino event rates obtained when overlapping the contribution from the different nuclear masses in the composite target.
We now turn to estimate the possible impact of NC-NSI on the predicted sensitivity of some proposed experiments. We will consider here four facilities: three with monoatomic detector targets (DARWIN, ARGO and Super-CDMS HV) and one using a molecular target (CRESST phase III). Their targets, energy thresholds and exposures can be found in Table 1. The results are shown in Fig. 4 where we plot our estimate of the range of variation of the neutrino floor induced by NC-NSI superimposed to the attainable bounds on the WIMP-nucleon spin independent cross section as reported by each of the experiments. It is important to stress that this figure cannot be used to make a final quantitative statement as the shown regions are calculated under different assumptions. On the one hand we compute the ultimate neutrino floor, i.e., considering only the neutrino irreducible background and for an experiment 100% efficient above its threshold. On the other hand the attainable bounds quoted by the experiments assume 30% (DARWIN), 90% (ARGO), 85% (Super-CDMS HV) and 70% (CRESST phase III) efficiency above their energy threshold. The curves also do not always correspond to the same CL. Nevertheless, the comparison can be used to foresee which experiments could potentially suffer a reduction of their sensitivity reach because of the presence of NC-NSI. As can be seen from this figure, for DARWIN [7] (top-left panel) the SM discovery limit is below the expected sensitivity, however after the inclusion of presently allowed NC-NSI the neutrino floor can rise above the experimental sensitivity for light WIMP masses m χ 20 GeV. This is also the case for the ARGO [9] (middle-left panel) and CRESST phase III [6] (bottom panel) experiments for which we find that NC-NSI can increase the neutrino floor above the expected sensitivity reach for WIMP masses m χ 50 GeV and m χ ∼ 4 − 10 GeV, respectively. Thus if these experiments detect a signal compatible with a WIMP around those mass ranges, it could be also interpreted as a neutrino event produced by flavor-dependent NC-NSI. Conversely the results show that these experiments have the potential to place meaningful bounds on NC-NSI. Note, however, that in all the cases presented here the floor has been somehow underestimated as a lower (more realistic) efficiency would move the discovery limits up, hence closer to the experimental reach.

Conclusions
Next generation direct detection DM experiments are expected to substantially increase the sensitivity to WIMP-nucleon cross sections for a wide range of DM masses. This will place them closer to the region where solar and atmospheric neutrinos undergoing CνNS will become an irreducible background.   [7], SuperCDMS -HV (Ge -top-right, Simedium-right panels) [8], ARGO (medium-left panel) [9] and CRESST phase III (bottom panel) [6]. The red dashed line corresponds to the SM prediction while the blue full line to the best fit NSI and oscillation parameters of the analysis of oscillation+COHERENT data. The shaded regions bordered by dashed, dot-dashed, and dotted lines correspond to predictions with NSI and oscillation parameters within ∆χ 2 OSC+COH ≤ 1, 4, 9, respectively. In each panel we show in green, for comparison, the sensitivity region as quoted by each of these experiments in the given references. For reference we also show by a green dashed line the region excluded by Xenon1T at 90% CL [87].
While the first measurement of CνNS by the COHERENT experiment observed scattering rates that seem to be in agreement with the SM predictions, BSM contributions of NC-NSI with the target nucleus, at a level that could impact future DM detectors, cannot be excluded. This is true even if neutrino oscillation constraints are taken into account.
We investigated in this work to what extent the future experimental discovery limits of DM direct detection facilities could be compromised by this flavor-dependent BSM. We used the full likelihood information of the global analysis of neutrino oscillation experiments and COHERENT scattering data to consistently account for the correlated effects of NC-NSI in the neutrino propagation in matter as well as in its interaction with the detector.
To be specific, we estimated the maximum effect which these still allowed BSM interactions could have in experiments with different targets (Xe in DARWIN, Ar in ARGO, Ge and Si in Super-CDMS HV, and CaWO 4 in CRESST phase III). We showed that the neutrino background contribution can be up to about five times larger than what is expected by the SM, certainly having impact on the ultimate reach of the ARGO, CRESST phase III and DARWIN experiments. In turn, these experiments may be able to constraint NC-NSI to yet lower levels than currently allowed by data.