Neutrino Backgrounds in Future Liquid Noble Element Dark Matter Direct Detection Experiments

Experiments that use liquid noble gasses as target materials, such as argon and xenon, play a significant role in direct detection searches for WIMP(-like) dark matter. As these experiments grow in size, they will soon encounter a new background to their dark matter discovery potential from neutrino scattering off nuclei and electrons in their targets. Therefore, a better understanding of this new source of background is crucial for future large-scale experiments such as ARGO and DARWIN. In this work, we study the impact of atmospheric neutrino flux uncertainties, electron recoil rejection efficiency, recoil energy sensitivity, and other related factors on the dark matter discovery reach. We also show that a significant improvement in sensitivity can potentially be obtained, at large exposures, by combining data from independent argon and xenon experiments.

The focus of this paper are direct detection experiments that look for local dark matter scattering on target nuclei in deep underground detectors [7][8][9]. Such experiments are particularly well-suited to finding weakly interacting massive particles (WIMPs) and WIMP-like dark matter [20][21][22][23][24]. These are among the most promising and best-studied candidates for dark matter. Most notably, they can obtain the observed dark matter relic abundance during the evolution of the early universe through the mechanism of thermal freeze-out [21]. New particles with masses near the weak scale are also predicted by most theories that address the electroweak hierarchy problem, and indeed many of them contain viable WIMP(-like) dark matter candidates [25,26].
Current direct detection bounds on WIMP dark matter candidates with masses above a few GeV (with mainly spin-independent interactions with nuclei) are dominated by experiments using liquid noble gases as their target material, primarily xenon and argon. Existing experiments have already achieved exposures of multi-tonne years while controlling backgrounds, as demonstrated in Refs. [27][28][29][30][31][32][33][34]. In coming years, the LZ [35], PandaX-4T [36], and XENONnT [37] are expected to achieve exposures near 20 t yr in xenon, while Dark-Side20k [38] aims for a total exposure near 100 t yr in argon. Beyond this, work is already under way on DARWIN based on xenon with an exposure goal of nearly 200 t yr [39], and ARGO based on argon reaching an exposure of 3000 t yr [40].
As such dark matter detectors grow in size and sensitivity, neutrinos from the sun [41,42], cosmic rays in the atmosphere [43][44][45], and diffuse supernovae [46][47][48] will become significant backgrounds to dark matter searches [49][50][51][52][53][54][55][56]. These neutrinos can scatter coherently with the target nucleus as a whole [57], as recently measured for the first time (for reactor neutrinos) in CsI(Na) [58] and argon [59]. In liquid noble element detectors this produces a nearly irreducible background to dark matter scattering on nuclei. Neutrinos can also scatter with electrons in the target material to produce a further background due to the finite ability of detectors to differentiate between nuclear and electronic interactions, especially close to the nuclear scattering detection energy threshold of a few keV [32,39].
If the neutrino fluxes and energy spectra were known precisely, neutrino scattering would present a challenging statistical background, reducing the scaling of sensitivity with total exposure to a square root rather than linear [52,53]. However, the full story is more complicated. On the one hand, differences between the recoil energy spectra from dark matter scattering relative to neutrinos, due to a combination of their respective particle physics cross sections and fluxes, provide a small degree of distinction for the signal. Unfortunately, on the other hand, uncertainties in the neutrino fluxes induce a more difficult systematic background that is hard to overcome when they grow larger than the size of the dark matter signal itself. These uncertainties can ultimately lead to a so-called neutrino floor beyond which it is difficult to make progress in the direct detection approach [52,53].
Future dark matter experiments based on xenon and argon, such as DARWIN [39] and ARGO [40], are expected to reach nearly down to their respective neutrino floors throughout almost the entire range of WIMP(-like) dark matter masses they are designed to look for.
Moreover, a broad range of WIMP(-like) dark matter theories predict nuclear scattering cross-sections near or below these detection-floors [60][61][62][63][64][65]. To optimize the dark matter reach of these experiments, it is therefore crucial to understand the neutrino backgrounds and predict how they depend on the properties of the detectors such as total exposure, electron rejection efficiency, and recoil energy sensitivity.
In this work, we investigate neutrino backgrounds to searches for WIMP(-like) dark matter through nuclear recoils in large-scale xenon and argon detectors. Our paper builds on previous work on neutrino backgrounds in three main ways [49][50][51][52][53][54][55][56]. First, we update detector energy sensitivity regions and resolutions for both future xenon and argon experiments, with a greater emphasis on argon relative to previous works. Second, we investigate the impact of electron recoil rejection on the sensitivity to dark matter subject to neutrinoelectron scattering backgrounds. And third, we show that combining data from argon and xenon experiments can improve the sensitivity to dark matter beyond just an increase in statistics. Along the way, we also study the effects of detector location, energy detection regions of interest, and spectral shape uncertainties on dark matter detection.
The outline of this paper is as follows. We begin in Section II by reviewing the calculations of the scattering rates for dark matter and background neutrinos in dark matter direct detection experiments. Next, in Section III we present the motivating assumptions we make about the detection capabilities of future detectors as well the statistical methods we use in our analysis. Section IV contains our main results on how the neutrino floor in xenon and argon is impacted by neutrino flux uncertainties, electron recoil rejection, and recoil energy sensitivity, as well as the effect of combining data from different future experiments. Our conclusions are presented in Section V. Some additional background material is summarized in Appendices A, B.

II. DARK MATTER AND NEUTRINO SCATTERING
In this section, we review the calculations of dark matter and neutrino scattering rates.
We also expand on previous treatments of neutrino scattering on electrons in detector targets.

A. Dark Matter Scattering
The rate R (N ) χ of dark matter χ scattering off a nuclear species N = (A, Z) per unit target mass is [9,23] d R where E R is the nuclear recoil energy, m χ is the dark matter mass, ρ χ is the local dark matter energy density, n N is the number of N nuclei per unit target mass, and dσ N /dE R is the differential scattering cross-section. The integral in Eq. (1) runs over the dark matter velocity v in the laboratory frame, with f lab ( v) being the lab-frame dark matter velocity distribution function. This integral is limited by v > v min , the minimum dark matter velocity needed to produce a recoil of energy E R . For elastic scattering, v min is given by with µ N = m χ m N /(m χ + m N ) the dark matter-nucleus reduced mass. Note as well that the number of N nuclei per unit detector mass n N can be expressed as where m N is the mass of a single N nucleus and F N is the mass fraction of N in the detector.
In this work, we focus on the broad range of WIMP and WIMP-like dark matter theories in which the dominant nuclear scattering interaction is spin-independent and mediated by massive intermediate particles. In this case, the differential cross-section takes the form [9,23] Hereσ N depends on the nuclear target but is independent of v and E R , and F N (E R ) is a nuclear form factor that is approximated well by [66] for s = 0.9 fm, and withā = 0.52 fm and c = (1.23 A 1/3 − 0.6) fm. To compare the sensitivities of dark matter detectors with different targets, it is also convenient to define a per-nucleon, spin-independent cross-section: where µ n = m n m χ /(m n + m χ ). Given the dark matter cross-section form of Eq. (4), the scattering rate of Eq. (1) can be written in the form The last term represents the integral over the dark matter halo and is given by We evaluate the halo integral using the Standard Halo Model (SHM) with the parameters v 0 = 238 km/s, v E = 254 km/s, and v esc = 544 km/s, as well as ρ χ = 0.3 GeV/cm 3 . These are the recommended values from Ref. [67] and they are in line with the values used in setting limits by recent direct dark matter search experiments [28][29][30]. In evaluating Eq. (9), we sum this expression over all isotopes of the target material weighted by their natural abundances.
Dark matter can also scatter off atomic electrons in the target material of the detectors.
However, DM-electron scattering produces a much weaker signal than the nuclear recoil searches we focus on in this study for two reasons. First, the energy transfer from weakscale dark matter to electrons is much less efficient when compared to scattering off nuclei, and very few electron scattering events produce visible energies above the nuclear recoil equivalent thresholds of E R 5 keV we focus on here [68][69][70]. Second, the searches we investigate reject electron recoils to a very high degree as a way to reduce backgrounds. For both reasons, we neglect dark matter-electron recoils in our analysis.

B. Neutrino Scattering on Nuclei
Coherent elastic scattering of neutrinos on nuclei was observed recently for the first time in CsI(Na) [58] and Ar [59] detectors, with rates that match Standard Model predictions. The same scattering processes from solar, supernova, and atmospheric neutrinos in dark matter detectors will become a significant background for future experiments as they become more sensitive. At leading order in the Standard Model, the differential cross-section for this process on a target nucleus N = (A, Z) comes from Z 0 vector exchange and is given by [71] where E ν is the incident neutrino energy, G F is the Fermi constant, F N (E R ) is the Helm form factor defined in Eq. (5), and with s 2 W = sin 2 θ W 0.23 [72]. Since this process is mediated by Z 0 boson exchange, it is independent of neutrino flavour. It is only allowed kinematically for E ν > E min with Therefore, given the nuclear scattering cross-section, the rate per unit target mass R where n N is the number of N nuclei per unit target mass, and the sum runs over all energydifferential neutrino fluxes Φ j (E ν ) with j = 1, . . . , n ν indexing the relevant sources. For each source, we have also implicitly summed over all neutrino flavours: for a = e, µ, τ,ē,μ,τ . In searches looking for heavier dark matter that induce nuclear recoils with E R ∼ 10-100 keV, the dominant nuclear recoil background comes from atmospheric neutrinos with energies in the range E ν ∼ 30-1000 MeV. Further details on neutrino sources relevant for dark matter searches are discussed in Sec. II D below.

C. Neutrino Scattering on Electrons
Neutrinos can also scatter elastically off of electrons. Such recoils in dark matter detectors can be misidentified as nuclear recoils. The probability for this to occur depends on the particle identification power of the experiment. Thus, these neutrino-induced electronic interactions can contribute to the background budget of dark matter searches based on nuclear scattering. The neutrino-electron scattering rate per unit target mass of neutral atoms containing the nucleus N = (A, Z) from neutrino species a = e, µ, τ,ē,μ,τ is where T is the energy transferred by the neutrino to the target and dσ (Ze) a /dT is the total differential cross-section for electron scattering on the atomic target.
To discuss the cross-section dσ (Ze) a /dT , it is convenient to begin with the cross-section for elastic neutrino scattering off a free electron, ν a + e − → ν a + e − , given by [73][74][75][76] dσ where the parameters Q + and Q − are and the minimum energy transfer E min is In many analyses of neutrino backgrounds for dark matter detection, neutrino scattering off electrons is computed using the free-electron approximation (FE) in which the atomic electrons of the target are treated as being unbound and at rest. The total neutrino crosssection off atoms of the nuclear species N = (A, Z) in the FE approximation is The intuitive motivation underlying this approximation is that the energy transferred by the neutrino to the electron is relatively large compared to (most of) the atomic binding energies, and thus atomic effects should not be important. Indeed, for neutrino-induced electron recoils with T = 1-100 keV relevant for weak-scale dark matter searches, the total scattering rate is dominated by solar neutrinos with energies between about E ν ∼ 0.1-1.0 MeV, with the largest contributions from the pp continuum and 7 Be line sources.
Going beyond the free-electron approximation, atomic effects on neutrino-electron scattering were studied in Refs. [77][78][79][80][81][82][83][84]. In Ref. [77], it was argued that atomic effects are small when E ν T, E n , where E n α 2 Z 2 m e /2n 2 are the relevant atomic binding energies. 1 This work was expanded in Refs. [78][79][80], which proposed a simple stepping approximation to capture the leading atomic effects in the regime E ν T, E n . This approximation extends the free-electron approximation formula of Eq. (20) by replacing Z with Z ef f (T ), given by where the sum runs over all relevant atomic levels n with binding energies E n . Careful ab initio calculations of the neutrino-electron differential cross-section in atoms were performed in Refs. [81,82] for germanium and in Ref. [83] for xenon. When considering pp and 7 Be solar neutrinos scattering on atomic electrons in xenon, Ref. [83] found that the full result is similar to the stepping approximation but smaller by about 25% up to T ∼ 10 keV and then falls off more quickly at higher energies. In the analysis to follow, we apply the results of Ref. [83] for neutrino-electron scattering on xenon but use the stepping approximation for neutrino-electron scattering on argon where it is expected to be a better approximation due to its lower binding energies (E n 3.2 keV in Ar versus E n 35 keV in Xe). We also assume that the energy transfer T produces visible energy in the detector with the same efficiency as free-electron recoils.
Atm (sum) Figure 1. Neutrino flux spectra relevant for dark matter direct detection with overall normalizations and source models as listed in Table I.

D. Neutrino Fluxes and Scattering Rates
The scattering kinematics discussed above imply that the neutrino backgrounds relevant for the direct detection of weak-scale dark matter correspond to neutrino energies E ν ∼ 1-1000 MeV, for neutrino-nucleus scattering and E ν ∼ 0.1-1 MeV for neutrino-electron scattering. In these energy ranges, the most important neutrino flux contributions come from atmospheric neutrinos created by cosmic ray showers [43][44][45]86], solar neutrinos from nuclear reactions in the Sun [41,42,[87][88][89], and the diffuse supernova neutrino background (DSNB) [46][47][48]. The neutrino flux spectra for the various sources are shown in Fig. 1 and are based on the source models listed in Table I. Each independent flux source j can be written in the form where φ j is the total flux and f j (E ν ) is an energy distribution normalized to unity. For reference, we list in Table I  include Refs. [43-45, 86, 90-93]. These calculations agree reasonably well with each other and with direct measurements of atmospheric neutrino fluxes by the Super-Kamiokande experiment down to energies of a few hundred MeV [94]. The estimated uncertainties on these calculations over the relevant energy range are roughly 20% [86], but with much smaller fractional uncertainties of about 5% in the flux ratios ν e /ν µ andν/ν [44,86,91]. In this work, we use the neutrino fluxes calculated in Ref. [44] for the simple reason that they are the most recent results that extend to the lower energies needed for this analysis, down to scenario of Ref. [97,98] as updated in Ref. [99]. This model appears to give a better agreement with observations than other proposals [96,100]. For the 8 B flux, we use the determination of Ref. [101] based on data from SNO [102] and SuperK [103] combined with fits to neutrino oscillation data. The 7 Be (862 keV) line flux adopted in this work comes from the recent measurement by Borexino [104]. We use the 7 Be branching ratios from Ref. [105] to fix the flux of the related 7 Be (384 keV) line in terms of the 7 Be (862 keV) flux.
The rate of neutrino-electron scattering relevant for dark matter searches also depends on 2 Relative to LNGS, the atmospheric fluxes computed in Ref. [86] for SNOLAB are larger by about 30%.
the electron-neutrino fraction F e of the solar neutrino flux from pp and 7 Be (862 keV); for both channels, we use the predicted and approximately constant MSW LMA [106][107][108] value F e = 0.55 (1 ± 0.02) based on recent determinations of the neutrino oscillation parameters collected in Ref. [72].
The collected emission from many supernovae over the history of the cosmos forms a diffuse supernova neutrino background (DSNB). This source is the dominant flux in the energy range E ν ∼ 20-40 MeV. To model the DSNB fluxes of the e,ē, and x = µ,μ, τ,τ flavours (at production), we apply the methods of Ref. [48]. In doing so, we follow Ref. [109] and use the parametrization of the star formation rate from Ref. [110] and connect it to the cosmic supernova rate based on Ref. [111]. The neutrino fluxes per supernova are based on Ref. [48] and use a simple Fermi-Dirac distribution with E tot ν = 3 × 10 53 erg, and effective temperatures T a = 4, 5, 8 MeV for a = e,ē, x. These estimates have a number of uncertainties associated with them, and in our analysis we apply an overall 50% uncertainty on the summed DSNB rate.
In Fig. 2 we show the total differential scattering rates per unit target mass as a function of recoil energy E R for neutrino-nucleus scattering in xenon (left) and argon (right). We also show the corresponding rates for dark matter scattering with cross-section σ n = 10 −48 cm 2 and masses m χ = 30, 100, 300 GeV. For xenon, we also show the rate for neutrino-electron scattering in terms of the reconstructed nuclear-recoil equivalent energy and scaled by a factor of 10 −3 , which is on the order of the expected electron rejection power of future xenon detectors. No such line is shown for argon since future detectors are expected to be able to reject electron recoils to a very high degree. The dashed, dotted, and dotdashed lines indicate the atmospheric (Atm), supernova (DSNB), and solar neutrino flux contributions to the neutrino-nucleus scattering rate. We see that the atmospheric flux is the dominant neutrino-nucleus background for E R 5 keV in xenon and E R 20 keV in argon. The recoil energy spectrum from atmospheric neutrinos is also similar to that from heavier m χ 100 GeV dark matter in xenon, while it is somewhat more distinct in argon.
Neutrino-electron scattering can become important in xenon at higher recoil energies, even with strong electron recoil rejection, but its spectral shape is much flatter than for dark matter at these energies.

III. DETECTING AND CHARACTERIZING SCATTERING SIGNALS
The signals induced by dark matter or neutrino scattering are shaped by the detectors used to measure them. In this section, we discuss the expected properties of future liquid noble element-based detectors and specify the assumptions we make to analyze their sensitivity. We also describe the statistical procedure used to test the discovery power of these experimental approaches.

A. Properties of Liquid Noble Element Detectors
Planned liquid noble element detectors will search for WIMP-like dark matter with mass above m χ 10 GeV primarily through nuclear recoils (NR). Many of the backgrounds to these searches come from electron recoils (ER). As such, these detectors are designed to have large acceptances for nuclear recoils together with the ability to reject the vast majority of electron recoils.
Future large-scale xenon detectors are expected to have a dual-phase format similar to that used in ZEPLIN [27], XENON [30,113,114], LUX [28,115], and Panda-X [29, 33,  116]. Scattering in the liquid bulk of the detector leads to scintillation, ionization, and heat. The primary scintillation is detected directly, while the ionization electrons are drifted through an electric field to a cathode at the liquid-gas interface region where they create secondary photons. The energy deposited by scattering events can thus be characterized by the number of primary photons S1 (scintillation) together with the number of secondary photons S2 (ionization). Combining S1 and S2, with appropriate weights after correcting for the position, gives the total energy deposited by the event. Finally, the ratio S2/S1 is used to discriminate between nuclear and electron recoils.
By extrapolating from the current capabilities of two-phase xenon detectors, future largescale detectors are expected to have excellent sensitivity to nuclear recoils together with strong electron rejection power for recoil energies above about E R 5 keV. Radioactive backgrounds to xenon dark matter searches tend to rise at higher recoil energies; to avoid the worst of them, the search region for dark matter scattering is typically limited from above to E R 35 keV. In the following analysis, we assume a region of interest (ROI) for weak-scale dark matter searches in xenon experiments of E R ∈ [5,35] keV, in line with the projections of Refs. [35][36][37]39]. Within this ROI, cutting on the S2/S1 ratio is expected to allow the rejection of electron recoils by a factor of 10 −4 -10 −3 at the cost of reducing the acceptance for nuclear recoils to about 30-50% [39]. For our analysis, we assume constant electron rejection factors over the entire ROI energy range between ε e = 2 × 10 −4 -2 × 10 −3 together with a realistic but optimistic nuclear recoil acceptance of ε n = 0.5. These properties are summarized in Table II.
Planned large-scale argon detectors [38,40] may rely on either a single-phase technique like DEAP-3600 [32,117] or a dual-phase format such as DarkSide [31,38]. The scintillation yield of primary photons in argon is smaller than in xenon (Xe ∼ 60 γ/ keV and Ar ∼ 40 γ/ keV [118,119]), and this tends to translate into higher energy thresholds for NR when relying on primary scintillation. Based on the recent performance of DEAP-3600 [32], we assume a nuclear recoil energy ROI in argon of E R ∈ [55,100] keV, although we also study the impact of reducing the lower energy threshold. A crucial feature of argon detectors is the ability to use pulse-shape discrimination to distinguish between nuclear and electronic interactions [120][121][122]. Based on the properties of current argon dark matter detectors, we consider flat electron rejection factors in the range ε e = 10 −9 -10 −7 along with a nuclear recoil acceptance over the ROI of ε n = 0.9 [32,38,123]. Similar to xenon, the adopted argon detector parameters are reported in Table II.
In both xenon and argon, electron recoils create more visible energy (and less energy lost to heat) than nuclear recoils for a given recoil energy [124]. As a result, an electron recoil that is misidentified as a nuclear recoil with energy E R corresponds to a smaller electron recoil energy T . The average relation between these quantities is defined as the quenching factor q ef f : This factor has been estimated using phenomenological models in Refs. [125][126][127][128], and investigated in direct calibration measurements in both argon [129][130][131][132] and xenon [133][134][135].
A selection of these experimental results have been used by the NEST collaboration to make a global fit for q ef f (E R ) [119,136]. In this work, we follow NEST and others and adopt a data-motivated model with the form where A q and B q are dimensionless coefficients. These parameters can depend on the strength of an applied electric field, which is relevant for dual-phase detectors. For xenon, the NEST collaboration quotes A q = 0.151 ± 0.027 and B q = 0.1 ± 0.05 at zero applied electric field [136] while the recent LUX analysis finds A q = 0.173 and B q = 0.05 over a range of applied fields between about 30-600 V/cm [137]. In the xenon analysis to follow we fix A q = 0.16 and B q = 0.08, assuming an E-field strength between 50-200 V/cm. For argon, in zero applied field, we set A q = 0.185 and B q = 0.12 which is consistent with the NEST estimate of A q = 0.19 ± 0.01 (for W = 15.3 eV) and B q = 0.101 ± 0.025 [136] as well as recent measurements by SCENE [130] and ARIS [131].
Realistic detectors also have a finite recoil energy resolution. In modern xenon and argon detectors, most of the spread σ R in the reconstructed recoil energy, E R , comes from stochastic fluctuations in the amounts of light and charge produced in individual scattering events [128,138]. To account for this factor in future xenon detectors, we assume a resolution similar to that achieved in the LUX [139,140] and XENON1T [141] experiments based on combined S1 and S2 energy determination. These experiments report an energy resolution for higher-energy electron recoils of σ/T 0.32/ T / keV. We extrapolate this resolution to nuclear recoils by including the quenching factor discussed above, The recent NEST determination of Ref. [136] argues that at electron energies below T ∼ 5 keV, the resolution σ/E 0.14 saturates at a constant value. However, other interpretations exist for the resolution of xenon at these lower energies, and we apply the form of Eq. (25) at all energies, consistent with Ref. [140]. For our projections in argon, we assume a nuclear recoil energy resolution equal to that achieved by the DEAP-3600 experiment [32], which over the energy ROI is fit accurately by the following parameterization: To account for detector energy resolution, we use a binned statistical analysis (described in detail below) with recoil energy bins constrained to be larger than the local energy resolution.
With an energy resolution of the general form σ R /E R = A r (E R / keV) −Br , as we find for both argon and xenon, the maximum number of bins satisfying this constraint is where x denotes the floor function and For i = 1, 2, . . . N b , the i-th recoil energy bin then covers the range This approach is expected to give a good approximation of the true energy resolution as long as the signal and background vary reasonably slowly from bin to bin, which we find to be the case in the examples studied here.
In addition to events from dark matter and neutrinos, experiments also detect background events from various sources. Projections for future experiments suggest that these backgrounds can be reduced to be much less important than neutrino-induced backgrounds over the dark matter search ROIs in argon and xenon detectors. In both detector materials, the dominant nuclear recoil background comes from neutrons produced by radioactive decays; these can be mitigated effectively by the use of radiopure detector materials and active veto systems [38-40, 138, 142]. Electron recoil backgrounds arise primarily from nuclear decays of noble element contaminants in the detector material. They are expected to be subleading relative to electron recoils from solar neutrinos in xenon [138], while they are reduced to low levels by purification and electron recoil rejection in argon [40]. Based on these considerations, we only consider backgrounds from neutrinos in our analysis.
Putting these detector considerations together, the mean detection rate per unit target mass for events identified as nuclear recoils with energy E R is All detector parameters relevant for evaluating this expression are collected in Table. II. For each bin i = 1, 2, . . . N b , the expected number N i of events is then where M T is exposure, equal to the total detector mass times the observation time.

B. Statistical Methods
Following Refs. [52,53,143,144], we use the profile likelihood method [145] to compute the dark matter discovery limit for a given set of experimental parameters, corresponding to the smallest dark matter cross-section for which the experiment can exclude the backgroundonly hypothesis with a 3σ significance at least 90% of the time. The test statistic for our analysis is a ratio of likelihoods that compares the background-only hypothesis to the signal hypothesis. As discussed above, we assume binned data organized in recoil energy bins i = 1, 2, . . . , N b . The expected number of events in bin i can thus be written in the form where ξ = σ/σ 0 is the χ-nucleon cross-section relative to some reference value and θ are a set of nuisance parameters that characterize the neutrino backgrounds. For a given set of data {N i }, the likelihood function is taken to be where P(n, λ) = λ n e −λ /n! is the Poisson distribution and L b ( θ) is a likelihood function for the background (and possibly signal) parameters. The test statistic for our analysis is then Here,ξ andˆ θ are the parameters that maximize the unconstrained likelihood of Eq. (33) for the set of data {N i }, and θ are the values that maximize it for the background-only hypothesis. Larger values of q 0 correspond to poorer fits for the background hypothesis to the data.
It is now straightforward to formulate the dark matter discovery limit in a precise way. We follow Refs. [52,53,143,144] and say that an experiment can discover a given dark matter scenario with cross-section strength ξ if it is able to exclude the background hypothesis with p-value below p < p 0 = 0.0013 (corresponding to an exclusion of at least 3σ for a Gaussian distribution) at least a fraction F = 0.90 of the time. Let us define q lim 0 by the relation The condition to exclude the background hypothesis is then To evaluate the dark matter discovery condition of Eq. (36) for a given scenario, we make use of a key feature of the test statistic of Eq. (34): its distribution has a simple asymptotic form for N tot = i N i 1. As discussed further below, this condition is found to hold for the scenarios of interest in this work for exposures M T 10 t yr. By applying Wilks' theorem [146,147], it was shown in Ref. [145] that the probability distribution for q 0 when the underlying data has value ξ is with and q 0,A being the value of the test statistic evaluated on the Asimov data set with N i = N i .
The dark matter discovery limit is then the smallest value of ξ for which q 0,A > 18.34.
In the analysis to follow, we apply this asymptotic relation to estimate how the dark matter discovery limit varies with the size and properties of hypothetical future experiments.
The approach requires only a single minimization over the space of θ per scenario. We do this minimization using a global pre-conditioning step followed by a local minimization based on the MMA algorithm [148] as implemented in the NLopt optimization package [149]. In contrast, a full treatment that does not use the asymptotic forms would require generating multiple sets of pseudo-data, each needing two minimizations, to estimate the distributions f (q 0 |ξ) and f (q 0 |0). Applied to simple examples, Refs. [145,[150][151][152][153] find that convergence is reasonably good for N tot 10. This condition is met for the main scenarios of interest in this work with exposure M T 10 t yr. We have also verified for several specific scenarios that the discovery limits derived with the asymptotic results based on the Asimov data set agree very well with the limits obtained with simulated pseudo-datasets.

C. Characterizing Uncertainties
To compute the likelihoods needed to obtain dark matter discovery limits, it is necessary to specify the background-parameter likelihood function L b ( θ) in Eq. (33). In this work, we focus on neutrinos as the primary background to dark matter. Other background sources can be mitigated such that they are much less important over the dark matter regions of interest that we study [38,39,138]. The main background uncertainties are therefore the neutrino fluxes from the n ν sources considered. For these, we mostly fix the shapes of the flux energy spectra f j (E ν ) and allow for variations only in the overall flux normalizations φ j , as defined in Eq. (22). It is therefore convenient to describe flux uncertainties with the variables with j = 1, . . . n ν andφ j being the central flux values listed in Table I. In terms of these variables, we take the background likelihood function to be a product of independent Gaussians for each source, where ∆θ j coincide with the fractional flux uncertainties listed in Table I.

IV. DARK MATTER DETECTION AND THE NEUTRINO FLOOR
We turn now to applying the methods discussed above to estimate the effect of neutrino backgrounds on the dark matter discovery sensitivities of future large-scale liquid argon and experiments. In particular, we investigate the impact of atmospheric flux uncertainties, electron-recoil rejection factors, and recoil energy sensitivity. We also study the potential gain in dark matter sensitivity that might be obtained by combining data from argon and xenon experiments.
For the baseline flux and detector parameters used in our analysis listed in Tables I and   II, the SI dark matter discovery limits we find are summarized in Fig. 3 Fig. 2). In the rest of the section, we investigate the effects of varying away from these baselines.

A. Impact of Atmospheric Fluxes and Uncertainties
Atmospheric neutrinos are the dominant source of neutrino-nucleus scattering background for weak-scale dark matter discovery. The nuclear recoil energy spectrum they induce can be very similar to dark matter scattering, particularly for certain specific dark matter masses.
The ability of future dark matter detectors to distinguish dark matter from atmospheric neutrinos through spectral shape information is therefore very sensitive to how well the atmospheric neutrino flux is known. Uncertainty on the net flux in the relevant energy range, summed over all flavours of neutrinos and antineutrinos, is estimated in Refs. [44,45,93] to be approximately 20%.
To demonstrate the impact of this uncertainty, we show in the upper panels of Fig The importance of the uncertainty in the atmospheric neutrino flux has a significant impact on the dark matter discovery reach in both xenon and argon detectors. The effect is important both for smaller dark matter cross-sections where the total event rate is dominated by atmospheric neutrinos, as well as for larger dark matter masses where the recoil energy spectrum from dark matter becomes similar to that from atmospheric neutrinos.
Improved measurements of the atmospheric neutrino flux at energies below E ν ∼ 100 MeV could therefore increase the reach of future large-scale dark matter searches. Projections for DUNE [154] suggest that a fractional uncertainty on the order of 10% may be achievable in this energy range. However, for DARWIN-or ARGO-scale detectors, the impact of this improvement appears to be fairly modest. An additional study of the effect of uncertainties in the shape of the atmospheric neutrino energy spectrum is presented in Appendix B with a similar conclusion.
Plotting the dark matter discovery reach as a function of M T , as shown in Fig. 4, also illustrates the structure of the neutrino floor phenomenon. As the neutrino background becomes so large that its uncertainty grows larger than the dark matter signal, increasing the exposure M T further yields almost no improvement. This effect is more pronounced in xenon than in argon since the signal and background energy spectra tend to be more similar for xenon. As the exposure is increased even further, to values beyond currently foreseeable capabilities, we find that the background shape is determined so well from data that the dark matter sensitivity starts increasing again, as √ M T . In effect, there are enough neutrino events in this regime for the background to be determined directly from the data rather than the input uncertainty. This "softness" of the neutrino floor was observed in Ref. [52,53] and studied in Ref. [56]. We will investigate it further below in relation to other

B. Impact of Electron Recoil Rejection
Neutrino-electron scattering is a background to dark matter-nucleus scattering due to the finite ability of detectors to distinguish between nuclear and electron recoils. The energy transfer from neutrinos to electrons is more efficient than for nuclear recoils, and thus it also leads to a significant decrease in nuclear recoil acceptance [138]. In contrast, at much larger exposures, a reduced electron rejection power can lead to significant degradation of the dark matter reach and solidification of the neutrino floor. So far, we have treated the electron rejection factor ε e as constant as a function of recoil energy with no uncertainty. In practice, however, this quantity will have an energydependent uncertainty both from the characterization of a given detector [137,155], as well as the challenge in relating the energies of electron recoils to nuclear recoils [136]. To investigate the potential effects of not knowing ε e exactly, we model this uncertainty in two simple ways.
First, we allow for a 20% energy-independent uncertainty relative to the central value of ε e .
And second, we allow ε e to vary independently in each energy bin by 10% relative to the central value of ε e to approximate the effect of uncertainty in the energy dependence of the electron rejection. For both uncertainty models, we treat the relevant variations as nuisance parameters in our profile likelihood analysis following the same approach as for the neutrino flux uncertainties.
Results for these two electron recoil rejection uncertainty models are given in Fig. 7, second model leads to a much more rigid neutrino floor, even for dark matter masses where the dark matter and neutrino-electron recoil spectra are significantly different. When the shape of electron recoil background is known, it can be determined very well with enough data from the highest energies measured, where it dominates over nuclear scattering (as shown in Fig. 2). It can then be extrapolated down to the lower energies where the dark matter signal is expected to lie. However, such an extrapolation is impeded by energydependent uncertainties since the electron scattering background at lower energies is no longer completely fixed completely by its value at higher energies.
In future experiments, it would be helpful for the collaborations to publish data-driven estimates for the uncertainties in electron recoil rejection as a function of the reconstructed recoil energy E R . This result illustrates a more general point that recoil-energy-dependent uncertainties reduce the ability of the profile likelihood to learn the background distributions.

C. Dependence on the Recoil Energy ROI
Our analysis has assumed dark matter nuclear recoil energy search regions motivated by projections for future large-scale noble element detectors. For xenon, this range is 5,35] keV while for argon it is E R ∈ [55,100] keV. In both cases, these regions of interest (ROIs) are chosen to maximize dark matter acceptance while accounting for expected detector sensitivities and mitigating potential (non-neutrino) backgrounds. Even so, it also potentially useful to understand the effect of expanding these search regions.
In Fig. 8 we show the effect on the dark matter discovery sensitivity of increasing the ROIs in xenon from E R ∈ [5,35] keV to E R ∈ [5, 100] keV (left), and in argon from E R ∈ [55,100] keV to E R ∈ [20, 100] keV. For xenon, the lower bound is limited by the flux of solar 8 B neutrinos. However, the upper bound could be raised with improvements in both the electron recoil rejection and nuclear recoil acceptance [138]. Even if this can be achieved with no significant increase in other backgrounds, it does not improve the dark matter discovery reach very much. At these larger recoil energies, both the dark matter and neutrino-nucleus scattering rates are strongly reduced by the small values of the nuclear form factor in this range. In contrast, expanding the ROI in argon down to E R as low as 20 keV produces a large increase in the sensitivity to dark matter. This is to be expected for lower dark matter masses below m χ 100 GeV where the maximum recoil energy is not very large, but the increase also extends to much larger masses. At such larger masses, a significant portion of the nuclear recoil spectrum lies below E R = 55 keV, so expanding the ROI at the low end also yields more dark matter events. To achieve a lower energy threshold in a (single-phase) argon detector would require a strong increase in the overall photo-collection efficiency, pushing to the statistical limitations of the intrinsic scintillating properties of argon [122,156].

D. Combining Data from Argon and Xenon Experiments
Dark matter and neutrino recoil energy spectra are significantly different in argon relative to xenon, as shown in Fig. 2. Furthermore, the dominant uncertainties in the neutrino backgrounds are strongly correlated between different experiments. These features suggest that combining data from future argon and xenon experiments might enhance the sensitivity to dark matter beyond just increasing the total effective exposure. We investigate this possibility here. A similar study using other detector materials was done in Ref. [53].
The natural generalization of the discovery test statistic discussed in Sec. III is to update where M i is the number of events in the i-th bin of the argon experiment and N j is the number of events in the j-th bin of the xenon experiment. The key feature of this combined likelihood is that the same neutrino background parameters apply to both sets of data. This is certainly true for solar and DSNB neutrinos, and also for atmospheric neutrinos provided the two detectors are at the same location. We also expect it to be a reasonable approximation for the atmospheric fluxes for detectors at different locations since their dominant uncertainties, the underlying interaction model and the primary cosmic ray flux spectrum, are strongly correlated [44,45,86,93].
The estimated result for the dark matter discovery reach of combining data from independent xenon and argon experiments is presented in Fig. 9   The enhancement of the dark matter discovery reach from combining argon and xenon data is striking at very large exposures. This combined sensitivity is much greater than either experiment individually. With only a single detector material, to a significant extent, it is possible to make up for the lack of dark matter in the background-only hypothesis by adjusting the neutrino flux parameters. When argon and xenon data are combined, this compensation effect is much less effective due to the different recoil energy spectra in these two materials owing to argon being much lighter than xenon. For exposures in line with proposed next-generation experiments, the increase in sensitivity from combining argon and xenon data is more modest. The increased combined sensitivity is also likely to be moderated by other background sources and uncertainties related to the individual detectors. Even so, the power of combining data from lighter and heavier detector materials offers a potential brute force approach to testing dark matter beyond the (not so rigid) neutrino floor.

V. CONCLUSIONS
In this paper, we investigated the nature of neutrino backgrounds for future large-scale liquid argon and xenon dark matter experiments. We reviewed the calculation of scattering rates of dark matter and neutrino scattering on nuclei and electrons in the target material in Section II. Next, in Section III we connected these scattering rates to the properties of expected future detectors and the statistical methods used to study them. We then applied these methods in Section IV to investigate the effect of neutrino scattering on the dark matter discovery sensitivity of such detectors.
Our results confirm previous studies of neutrino backgrounds to dark matter discovery in xenon detectors based on a profile likelihood analysis [52][53][54]56]. As these backgrounds approach the dark matter signal being looked for, the scaling of the discovery sensitivity with exposure M T goes from linear to square root. When the systematic uncertainties in the background rates approach the size of the signal being looked for, the scaling with exposure slows even further. This is commonly referred to as the neutrino floor. However, if the recoil energy spectra of these backgrounds are known precisely, with enough additional data, a dark matter signal can be distinguished from neutrinos through spectral information and scaling with the square root of exposure resumes [52,53,56].
The most important neutrino background in searches for WIMP(-like) dark matter with mass m χ 10 GeV comes from atmospheric neutrinos. The search reach of future detectors is therefore sensitive to how well the total atmospheric flux can be determined, particularly in the neutrino energy range E ν ∈ [30, 300] MeV. We find that direct measurements of this flux could be extremely valuable to future dark matter searches. Updated calculations of the atmospheric neutrino flux as functions of the time and location of dark matter searches would also be helpful.
A further important source of backgrounds is sub-MeV energy solar neutrinos scattering off atomic electrons in dark matter detectors. Pulse shape information can be used in argon to distinguish this source from dark matter scattering on nuclei [32]. Distinguishing electron recoils is a greater challenge in xenon, and future detectors are expected to have electron rejection factors in the range ε e 2 × 10 −4 -2 × 10 −3 [35][36][37]39]. Our results suggest that this level of rejection is largely sufficient to handle neutrino-electron scattering backgrounds for the current and next generation of xenon detectors with total exposures up to M T ∼ 200 t yr. This is possible, in part, because the recoil energy spectrum from neutrino-electron scattering is significantly different from that of nuclear recoils, allowing for this background to be determined from data. However, we also note that if the shape of the background is uncertain, its determination from data becomes much more difficult.
Detectors can also be characterized by the recoil energy range over which they are sensitive to dark matter. For xenon, we assumed a region of interest for nuclear recoils of E R ∈ [5,35] keV. This range does a good job of covering dark matter masses above about 10 GeV; extending this region to lower energies introduces a large neutrino background from 8  As demonstrated in Eq. (9), all the dependence of the dark matter scattering rate on the astrophysical dark matter distribution can be collected into the local dark matter energy density ρ χ and the so-called halo integral [9,23], where v is the dark matter velocity in the lab frame, f lab is the lab-frame velocity distribution, and v min (E R ) = m N E R /2µ 2 N for a given recoil energy E R . The lab-frame velocity distribution is related to the local distribution in the galactic halo frame f ( v ) by where v is the dark matter velocity in the halo frame, and v E is the net velocity of the Earth relative to the dark matter halo.
To evaluate the halo integral, we follow the recent recommendations of Ref. [67] and set ρ χ = 0.3 GeV cm −3 for the local dark matter density together with the Standard Halo Model (SHM) velocity distribution, The factor of N is a normalization factor given by is the usual error function. In our analysis we use the recommended SHM parameters from Ref. [67]: v 0 = 238 km/s, v E = 254 km/s, and v esc = 544 km/s.
A convenient feature of the SHM is that it allows for an analytic expression for the halo integral defined in Eq. (10).

Appendix B: Shape Uncertainties in the Atmospheric Flux Spectrum
In the analysis above, we showed that atmospheric neutrinos are typically the most important contribution to the neutrino background to dark matter detection. This analysis assumed an uncertainty in the normalization of the atmospheric neutrino energy spectrum with the spectral shape held fixed based on the calculation of Ref. [44]. However, the energy dependence of the atmospheric flux also has an uncertainty. We investigate the effect of such an uncertainty on dark matter discovery in this appendix.
As a first step, it is helpful to connect specific atmospheric neutrino energies to the nuclear scattering recoil energies they induce. Collecting neutrino energies according to E k = 10 MeV 10 4 MeV 10 MeV we define k = 1, 2, . . . N ν = 20 neutrino energy bins with range E ν ∈ [E k−1 , E k ]. With these in hand, consider now the quantities For each recoil energy bin i, µ ik is the contribution to that recoil bin per unit exposure from atmospheric neutrinos with energies in the k-th neutrino bin. Multiplying by exposure M T would then produce the expected number of events from these energy ranges. In Fig. 10 we show the partial counts per unit exposure µ ik for xenon (left) and argon (right). The detector parameters used are the same as used previously and summarized in Table II. This figure indicates that the dominant contribution to nuclear scattering from atmospheric neutrinos comes from the neutrino energy range E ν ∼ 50-200 MeV. This is largely due to the flux peaking in this range, but also from the kinematics of the scattering covered in Sec. II B.
To estimate the impact of atmospheric spectral shape uncertainties on the dark matter discovery sensitivity of future detectors using the profile likelihood framework, we need a likelihood function for variations in the shape function. Recall that previously we only varied the overall atmospheric flux normalization with a Gaussian likelihood weight given in Eq. (41). If we had a motivated functional parametrization for the atmospheric spectral, it would be straightforward to generalize this likelihood function to include variations in the model parameters. However, since we do not have such a form, we pursue a different but motivated approach.
Our generalized likelihood function for the atmospheric flux shape is with N a constant normalization and S is given by where y = ln(E ν /10 MeV), y max = ln(10 GeV/10 MeV), Z is a constant, is the local fractional variation of the atmospheric neutrino flux spectrum relative to the fiducial value from Ref. [44], and ∆θ = 0.2 is the estimated fractional uncertainty on the total atmospheric flux. The first term in Eq. (B4) reproduces the previous likelihood when θ is constant in y. The second term is new and depends on the constant Z; it imposes a penalty on local variations in θ(y). For example, for Z = 1/4 a linear variation in θ(y) from θ(0) = −1 to θ(y max ) = +1 would contribute the same likelihood cost (from the derivative term) as the global variation θ(y) = 1 everywhere (from the non-derivative term). While the form of Eq. (B4) is somewhat arbitrary, we argue that, for Z ∼ 1, it allows for a reasonable estimate of the impact of atmospheric flux shape uncertainties on dark matter discovery. Ultimately, a direct comparison to atmospheric flux measurements such as in Ref. [94] would be preferable provided they can be extended to lower neutrino energies.
Our strategy to implement local variations in the atmospheric flux in estimating the dark matter discovery sensitivity is to use a discretized form of θ(y). As in Eq. (B1) above, we split the atmospheric neutrino energies into logarithmic bins k = 1, 2, . . . N ν = 20 ranging between E ν ∈ [10 1 , 10 4 ] MeV. The contribution to the number of atmospheric neutrino events in recoil energy bin i is where θ k is a (weighted) discretized form of θ y that characterizes the variation in the atmospheric flux over E ν ∈ [E k−1 , E k ]. These θ k variables are then treated as nuisance parameters in the profile likelihood with a likelihood weight motivated by Eq. (B4) of where we set Z = 1/4. To achieve a more robust minimization in the profile likelihood, we find it convenient to first minimize for rigid θ = θ 1 = θ 2 = . . . = θ Nν , and then reminimize in variations around this constrained minimum.
The results of this approach for the dark matter discovery limits in xenon and argon are shown in Fig. 11. Allowing for an increased freedom in the spectral shape of the atmospheric neutrino flux degrades the sensitivity to dark matter, but only by a very small amount. This indicates that allowing the normalization of this flux source to float captures most of the uncertainty on it. The result is also consistent with the doubly differential event rates shown in Fig. 10, where we see that the partial atmospheric event counts in each nuclear recoil energy bin are dominated by atmospheric neutrinos with energies close to E ν ∼ 100 MeV.