Impact of form factor uncertainties on interpretations of coherent elastic neutrino-nucleus scattering data

The standard model coherent elastic neutrino-nucleus scattering (CE$\nu$NS) cross section is subject to nuclear form factor uncertainties, mainly driven by the root-mean-square radius of the neutron density distribution. Motivated by COHERENT phases I-III and future multi-ton direct detection dark matter searches, we evaluate these uncertainties in cesium iodide, germanium, xenon and argon detectors. We find that the uncertainties become relevant for momentum transfers $q\gtrsim 20$ MeV, are essentially independent of the form factor parameterization and are larger for light elements, reaching up to 11% for argon. Consequently, form factor uncertainties are not important for CE$\nu$NS induced by reactor or solar neutrinos. Taking into account these uncertainties, we then evaluate their impact on measurements of CE$\nu$NS at COHERENT, the diffuse supernova background (DSNB) neutrinos and sub-GeV atmospheric neutrinos. We show that the number of events in COHERENT involves relative uncertainties that---depending on recoil energy and nuclide---can be as large as 52%. For DSNB and atmospheric neutrino fluxes the uncertainties can reach values up to 7% and 16%, respectively. Finally, we consider their impact on searches for nonstandard neutrino interactions, sterile neutrinos and neutrino generalized interactions. We point out that studies of new physics using CE$\nu$NS data are greatly affected by neutron form factor uncertainties, which if not properly taken into account may lead to the misidentification of new physics signals. The uncertainties quantified here are also relevant for dark matter direct detection searches.


I. INTRODUCTION
Coherent elastic neutrino-nucleus scattering (CEνNS) was observed by the COHERENT experiment in 2017 in a cesium iodide scintillation detector. The measurement used neutrinos produced at the spallation neutron source at the Oak Ridge National Laboratory [1]. The cross section, obtained by the coherent sum of the individual nucleon amplitudes, is the largest of all neutrino cross sections at energies E ν 100 MeV, exceeding the elastic neutrino-electron scattering cross section by about two orders of magnitude in typical nuclei. The observation, however, relies on the detection of very small recoil energies, which only recently became possible with the use of the technology employed in direct detection dark matter (DM) searches.
CEνNS data allow precise measurements of the weak mixing angle [2], detailed studies of nuclear structure through weak neutral current interactions [3] and opens the possibility of searching for new physics beyond the standard model (SM) [4,5]. Indeed, since its observation various studies of nonstandard neutrino interactions (NSI) [6][7][8][9], sterile neutrinos [9], neutrino generalized interactions (NGI) [10] and neutrino electromagnetic properties [8,11] have been presented. A proper interpretation of a CEνNS signal, as related to any of these new physics scenarios, requires a detailed understanding not only of experimental systematics errors but also of theoretical uncertainties.
The calculation of event rates in CEνNS experiments involves proton and neutron nuclear form factors, which account for the proton and neutron distributions within the nucleus. In most treatments, however, the form factors are assumed to be equal and so the nuclear form factor becomes a global factor, which is typically parametrized in terms of the Helm [12] form factor, the Fourier transform of the symmeterized Fermi distribution [13], or the Klein-Nystrand form factor [14] as adopted by the COHERENT collaboration [1]. These form factors depend on various parameters whose values are fixed via experimental data and so involve experimental uncertainties. Of particular relevance is the root-meansquare (rms) radius of the nucleon distribution, which in such analyses is fixed by, for example, the value derived using a particular nuclear physics model or through a value derived from fits to nuclear data [15]. This simplification introduces an uncertainty on the predicted CEνNS recoil spectrum (and number of events) in the SM as well as in beyond-thestandard-model (BSM) physics scenarios.
The root-mean-square (rms) radius of the proton distribution r p rms is known from elastic electron-nucleus scattering with a precision of order one-per-mille for nuclear isotopes up to Z = 96 [16]. This is in sharp contrast with the rms radius of the neutron distribution r n rms , which for almost all nuclear isotopes is poorly known. Theoretical uncertainties on the CEνNS process are therefore driven by the uncertainties in r n rms . For the COHERENT experiment, the quenching factor and neutrino flux uncertainties are of order ∼40% [1,17]. Thus, form factor uncertainties are not particularly relevant in the interpretation of current data. This situation, however, is expected to change in the near future, and so form factor uncertainties will play an important role.
Identifying the size of these uncertainties is crucial for two reasons: (i) To understand whether a given signal is the result of new physics or of an "unexpected" nuclear physics effect, (ii) DM direct detection in multi-ton detectors like DarkSide-20k [18], will be subject to irreducible neutrino backgrounds from the diffuse supernova background (DSNB) and sub-GeV atmospheric neutrino fluxes. A precise understanding of this background is crucial to discriminate between neutrinoinduced and WIMP-induced signals.
With this in mind, in this paper we investigate the size and behavior of the neutron nuclear form factor uncertainties and their impact on the interpretation of data. To that end we consider four well-motivated nuclides: cesium iodide, germanium, xenon and argon. The first three are (or will be) used by COHERENT in one of its three phases [19], while argon will be used by DarkSide-20k which will take 1000 ton-year of data [18]. For definitiveness we consider three nuclear form factor parametrizations: The Helm form factor [12], the Fourier transform of the symmeterized Fermi distribution [13] and the Klein-Nystrand form factor [14]. And we assume the same parametrization for both protons and neutrons. We first study the size and momentum transfer (q) dependence of the neutron form factor uncertainties using these three parametrizations. After precisely quantifying them, we study their impact in COHERENT and in an Argon-based multi-ton DM detector. We assess as well the impact of the uncertainties on the interpretation of new physics effects. We do this in the case of NSI, active-sterile neutrino oscillations in the 3+1 framework and spin-independent NGI. We evaluate the effects of the neutron form factor uncertainties on the available parameter space and the potential misidentification of new physics signals, when these uncertainties are not properly accounted for.
The paper is organized as follows. In Section II we introduce our notation and briefly discuss the CEνNS process, focusing on form factor parameterizations and the corresponding rms radii of the nucleon density distributions. In Section III we quantify the size of the neutron form factor uncertainties, study their q dependence and show that they are fairly independent of the choice of the nuclear form factor. In Sections IV and V we study the implications for SM predictions and for new physics searches, respectively. In Section VI we present our conclusions. In Appendix A we present the details of the calculation of the DSNB neutrino flux, while in Appendix B we provide details of the NGI analysis.

II. COHERENT ELASTIC NEUTRINO-NUCLEUS SCATTERING
For neutrino energies below ∼ 100 MeV the de Broglie wavelength of the neutrino-nucleus process is larger than the typical nuclear radius and so the individual nucleon amplitudes add coherently. In the SM this translates into a cross section that is approximately enhanced by the number of constituent neutrons N [20,21]: where E r = q 2 /2m N is the nuclear recoil energy. This result follows from the vector neutral current. The axial current contribution, being spin dependent, is much smaller. The neutron and proton charges are given by g n V = −1/2 and g p V = 1/2 − 2 sin 2 θ W , with θ W the weak mixing angle. In the Born approximation, the nuclear form factors F N,Z (q 2 ) follow from the Fourier transform of the neutron and proton distributions. They capture the behavior one expects: the cross section should fall with increasing neutrino energy (increasing q). Theoretical predictions based on Eq. (1) involve uncertainties from electroweak parameters and nuclear form factors. These uncertainties should be accounted for and are particularly important in searches for new physics effects, which arguably are not expected to significantly exceed the SM expectation. Since the uncertainty in G F is a few tenths of a part per million [22], electroweak uncertainties are dominated by the weak mixing angle for which (using the MS renormalization scheme at the Z boson mass scale) the PDG gives [23] Electroweak uncertainties are therefore of no relevance. On the contrary, since nuclear form factors encode information on the proton and neutron distributions one expects these uncertainties to be sizable and more pronounced for large q, given the behavior F(q 2 → 0) → 1. These uncertainties turn out to be crucial for the interpretation of data from fixed target experiments such as COHERENT [1,17] and for DM direct detection experiments subject to diffuse supernova background (DSNB) and sub-GeV atmospheric neutrinos [24].

A. Nuclear form factors
Form factors are introduced to account for the density distributions of nucleons inside the nucleus. They follow from the Fourier transform of the nucleon distributions, F(q 2 ) = e i q· r ρ(r) d 3 r . ( The basic properties of nucleonic distributions are captured by different parametrizations. Here we consider those provided by the Helm model [12], the symmeterized Fermi distribution [13] and the Klein-Nystrand approach [14]. These distributions depend on two parameters which measure different nuclear properties and which are constrained by means of the rms radius of the distribution, In what follows we briefly discuss these parametrizations and the relations between their defining parameters and the rms radius of the distributions. These relations are key to our analysis for they determine, through r rms experimental uncertainties, the extent up to which these parameters can vary, thereby defining the form factor uncertainties. In the Helm model the nucleonic distribution is given by a convolution of a uniform density with radius R 0 (box or diffraction radius) and a Gaussian profile. The latter is characterized by the folding width s, which accounts for the surface thickness. Thus, the Helm distribution reads with θ(x) a Heaviside step function and f G (x) a Gaussian distribution given by The Helm form factor is then derived from Eqs. (3) and (5): where j 1 (x) = sin(x)/x 2 − cos(x)/x is a spherical Bessel function of order one. The rms radius is obtained from Eqs. (4) and (5) and is given by The symmeterized Fermi density distribution ρ SF (r) follows from the symmeterized Fermi function f SF (r) = f F (r) + f F (−r) − 1, which in turn follows from the conventional Fermi, or Woods-Saxon function, where c is the half-density radius and a represents the surface diffuseness. Accordingly, ρ SF (r) can be written as ρ SF (r) = 3 4πc(c 2 + π 2 a 2 ) sinh(c/a) cosh(r/a) + cosh(c/a) .
In contrast to the Fermi density distribution it has the advantage that its Fourier transform can be analytically evaluated with the result, Then, The Klein-Nystrand approach relies on a surface-diffuse distribution which results from folding a short-range Yukawa potential with range a k , over a hard sphere distribution with radius R A . The Yukawa potential and the hard sphere distribution can be written as The Klein-Nystrand form factor can then be calculated as the product of two individual Fourier transformations, one of the potential and another of the hard sphere distribution, resulting in The rms radius is given by

III. FORM FACTOR UNCERTAINTIES
The rms radii of the proton density distributions are determined from different experimental sources. The values reported in [16] include data from optical and K α X-ray isotope shifts as well as muonic spectra and electronic scattering experiments. This wealth of data has allowed the determination of r 2 p ≡ r p rms with high accuracy for all isotopes of interest for CEνNS and DM direct detection experiments. The rms radii for the proton distribution are as in Table I. In contrast, rms radii of the neutron density distributions r 2 n ≡ r n rms are poorly known, mainly because barring the cases of 208 Pb, 133 Cs and 127 I [25][26][27], their experimental values follow from hadronic experiments which are subject to large uncertainties. At the form factor level, therefore, uncertainties on r p rms are basically irrelevant while uncertainties in r n rms have a substantial effect. Consequently, we adopt the following procedure. We verified that adopting different form factors for the proton and neutron distributions leads to a small effect on our results, so we assume the same form factor for both. For protons, in each of Eqs. (8), (12) and (15), we fix one parameter and determine the other by fixing r p rms to its experimental central value. For neutrons we do the same as that for protons, but allow r n rms to takes values above r p rms ; this lower limit is reliable provided N > Z, which is the case for all the nuclides we consider. We fix the upper limit using the neutron skin, ∆r np = r n rms − r p rms , of 208 Pb, which is measured by the PREX experiment at Jefferson laboratory to be ∆r np ( 208 Pb) = 0.302 ± 0.177 fm [26]. Requiring that the neutron skin of the nuclides under consideration be no larger than 0.5 fm is reasonable given that they have N < N( 208 Pb). We use For the Helm form factor we fix the surface thickness s to 0.9 fm [15], for the form factor based on the symmeterized Fermi function we fix the surface diffuseness a to 0.52 fm [28] and for the Klein-Nystrand form factor we fix the range a K of the Yukawa potential to 0.7 fm [14]. We checked that our results are rather insensitive to variations of these values. With the procedure already outlined, we first investigate the behavior of the uncertainties and their size. Figure 1 shows the result for the Helm form factor obtained for 133 Cs. The left panel shows that for low q, the uncertainties are small and increase with increasing momentum, reaching their maximum  TABLE I. Rms radii (in fm) of the proton density distributions of the stable isotopes of cesium, iodine, argon, germanium and xenon [16]. Relative abundances of Ar, Ge and Xe are provided in parentheses.   Table I for  for q 65 MeV. This behavior is apparent in the middle panel which shows the Helm percentage uncertainty, which measures the size of the spread due to the uncertainties in r n rms . It can be seen that in the case of 133 Cs the uncertainty can be as large as 8%, and for 40 Ar as large as 11%. To address how this result depends on the choice of form factor, we calculate the percentage uncertainty for F SF and F KN , with the aid of Eqs. (12) and (15). The right panel in Fig. 1 shows the relative uncertainty obtained by comparing the uncertainties from the Helm and the symmeterized Fermi form factors, calculated according to |U SF − U H |/U SF × 100%; results using the Klein-Nystrand form factor are similar and are not displayed. It can be seen that uncertainties are parametrization independent for q up to 60 MeV or so. For larger q, differences are at most of order 2.5%, with the Helm form factor yielding slightly larger values. In summary, the conclusions derived from Fig. 1 hold no matter the form factor choice. Henceforth, to calculate the impact of the form factor uncer-tainties on CEνNS, we employ the Helm form factor.

IV. IMPLICATIONS FOR COHERENT, DSNB AND SUB-GEV ATMOSPHERIC NEUTRINOS
We now turn to the study of the impact of the form factor uncertainties on SM predictions for CEνNS. We begin with COHERENT in each of its phases. For phase-I we calculate the expected number of events taking into account the contributions from both 133 Cs and 127 I. For phase-II (germanium phase) and phase-III (LXe phase) we calculate the number of events assuming the specifications given in [19] with the number of protons on target n POT fixed as in the CsI case (event numbers for a different value n POT can be straightforwardly rederived by scaling our result by n POT /n POT ). Since germanium and xenon have several sufficiently abundant isotopes (see Table I), we calculate the recoil spectrum generated by each of the nuclides. The i th isotope recoil spectrum can be  Table I. The Helm form factor has been used but a different form factor does not noticeably alter the spectra.
written as where m det refers to the detector mass, X i to its relative abundance, m to the average molar mass calculated as ∑ k m k X k (m k being the molar mass of the individual isotopes), N A = 6.022 × 10 23 mol −1 and φ(E ν ) the neutrino flux. Note that the global factor (m det N A / m )X i corresponds to the number of nuclei of the i th type in the detector. The differential cross section is given by Eq. (1) with m N → m i , N → N i and q 2 → q 2 i = 2m i E r . For each isotope contribution r p rms is fixed according to the values in Table I and r n rms as described in the previous section. Calculating the individual recoil spectra according to Eq. (18) and then summating over all of them (to determine the total recoil spectrum), allows to properly trace the uncertainties induced by each neutron form factor.
For the COHERENT phase-I analysis we use m det = 14.6 kg and adapt Eq. (18) to take into account the contributions of 133 Cs and 127 I. This is done by trading X i for the nuclear fractions f i = A i /(A Cs + A I ) (A i refer to the 133 Cs and 127 I mass numbers) in (18) and m for m CsI = 259 × 10 −3 kg/mol (CsI molar mass). The acceptance function is [17] A(n PE ) = where k 1 = 0.6655, k 2 = 0.4942, x 0 = 10.8507, and n PE is the observed number of photoelectrons. 1 Neutrino fluxes in CO-HERENT are produced by π + and µ + decays, and so three neutrino flavors are produced (ν µ ,ν µ and ν e ) with known energy spectra: The neutrino flux per flavor φ a (E ν ) at the detector is obtained by weighting the energy spectra by the normalization factor N = r × n POT /4πL 2 . Here r = 0.08 determines the number of neutrinos produced per proton collision (per flavor), L = 19.3 m is the distance from the source to the detector, and n POT = 1.76 × 10 23 is the number of protons on target in the 308.1 days of neutrino production [1], which corresponds to 2.1 × 10 23 protons/year. In terms of n PE we calculate the SM expectation for the number of events per bin (2 photoelectrons) taking into account the neutron form factor uncertainties. The result is displayed in Fig. 2. Uncertainties in the neutron form factor produce an uncertainty in the expected number of events, with a behavior such that small values of r n rms tend to increase the number of events, while large values tend to decrease them. This is in agreement with the result in the left panel in Fig. 1. One can see as well that for low n PE (recoil energy), no sizable uncertainties are observed. However, for n PE = 7 (E r = 5.98 keV) uncertainties are already of order 7% and increase to 16% for n PE = 15. We also calculate the number of events by fixing the rms radii of the 133 Cs and 127 I neutron density distributions to r n rms ( 133 Cs) = 5.01 fm , r n rms ( 127 I) = 4.94 fm . (21) These values follow from theoretical calculations using the relativistic mean field (RMF) NLZ2 nuclear model [27]. The result obtained can then be regarded as a purely theoretical result. Comparing the black-dotted histogram in Fig. 2 with those determined by the form factor uncertainties we see that the theoretical expectation is closer to the result for r p rms . COHERENT phase-II consists of a p-type point-contact high purity germanium detector with m det = 15 kg and located at L = 22 m from the source. COHERENT phase-III, instead, aims at measuring CEνNS by using a two-phase liquid xenon detector with m det = 100 kg and located at L = 29 m. At low recoil energies the number of CEνNS events in the Xe detector will exceed those in the Ge detector by about an order of magnitude. However, since Ge isotopes are lighter than Xe isotopes, the Ge detector will be sensitive to CEνNS events at higher recoil energies and so they are complementary [19]. Using these target material masses, the corresponding locations and assuming n POT as in the CsI calculation, we calculate the impact of the uncertainties on the expected number of events in both detectors. We adopt the same acceptance function as for the CsI detector. As can be seen from Fig. 3, in the germanium case, relative uncertainties are never below 2.5%, and increase up to 23% for large recoil energies. For Xe, relative uncertainties are larger with values as large as 52% for large recoil energies and never smaller than 5%. It is clear then that in the analysis of COHERENT data, form factor uncertainties should be taken into account.
Since form factor uncertainties increase with increasing momentum transfer, they are also relevant for CEνNS induced by the DSNB and sub-GeV atmospheric neutrinos. DSNB neutrinos (neutrinos and antineutrinos of all flavors) result from the cumulative emission from all past core-collapse supernovae. Their flux is thus determined by the rate for corecollapse supernova (determined in turn by the cosmic star formation history), and the neutrino emission per supernova, properly redshifted over cosmic time [29] (see Appendix A for details). The latter is well described by a Fermi-Dirac distribution with zero chemical potential and with T ν e < Tν e < T ν x [30]. For the calculation of the DSNB neutrino flux we use T ν e = 3 MeV, Tν e = 5 MeV and T ν x = 8 MeV, and sum over all flavors.
Atmospheric neutrino fluxes (ν e and ν µ and their antiparticles) result from hadronic showers induced by cosmic rays in the Earth's atmosphere. We take the atmospheric fluxes from Ref. [31] generated by a FLUKA Monte Carlo simulation [32], that includes ν e,τ andν e,τ fluxes up to about 10 3 MeV. We only consider atmospheric neutrino fluxes below 100 MeV because for higher energies the loss of coherence for CEνNS drastically depletes the neutrino event rate making the flux at those energies less relevant. Figure 4 shows the recoil spectrum obtained in an argon detector assuming an exposure of 5 × 10 3 ton-year. The dashed (solid) curve is obtained by fixing r n rms = r p rms (r n rms = 6 fm). In the calculation we include only 40 Ar and checked that form factor uncertainties for the "high energy" tail of the solar neutrino spectrum ( 8 B and hep neutrinos) are not relevant, as expected from the middle panel in Fig. 1. Up to recoil energies of order 18 keV the recoil spectrum is by far dominated by hep and 8 B neutrinos. DSNB-induced recoils thus dominate the spectrum in a small recoil energy window up to the recoil energy at which atmospheric neutrinos become important. The DSNB contribute sizably to the recoil spectrum up to energies of order 30 keV. At that energy we find that of the total number of events, the DSNB contributes about 18%. Note that form factor uncertainties for the recoil spectrum induced by DSNB neutrinos range from 4% at E r = 20 keV to 7% at E r = 30 keV. For atmospheric-induced recoils, instead, these uncertainties range from 4% for E r = 20 keV to 16% for E r = 70 keV, where about one event is expected.

V. IMPLICATIONS FOR NEW PHYSICS SEARCHES
We now discuss the effects of the neutron form factor uncertainties on the predictions for new physics. To do so, we consider three new physics scenarios that have been discussed in the literature in connection with COHERENT data: NSI [7,8], sterile neutrinos [8] and NGI [10].

A. Theoretical basics
NSI is a parametrization of a new physics neutral current interaction mediated by a vector boson of mass m V [33]. Dropping the axial coupling, which yields nuclear spin-suppressed effects, and in the effective limit (m V q CEvNS ), Written this way, the NSI parameters measure the strength of the new interaction compared to the weak interaction, ε q i j g 2 qi j /m 2 V /G F , where g qi j are gauge couplings. In the presence of NSI, the differential cross section becomes lepton-flavor dependent. For the i th neutrino flavor it can be derived from Eq. (1) by trading g n V → g n V + ε n i j and g p V → g p V + ε p i j , with ε n i j = ε u i j + 2ε d i j and ε p i j = 2ε u i j + ε d i j [4,5]. Oscillations with an eV mass sterile neutrino have an effect on the CEνNS event rate. If the flux of ν α neutrinos (α = e, µ, τ) at the source is Φ ν α (E ν ), the flux at the detector will be diminished by the fraction of neutrinos that oscillate into the sterile and the other active states. Quantitatively this means that the flux of neutrinos of flavor α that reach the detector is P αα Φ ν α (E ν ), where P αα is the ν α survival probability defined as P αα = 1−P αs −P αβ (α = β), with P αs and P αβ the ν α → ν s and ν α → ν β neutrino oscillation probabilities, respectively.
For short-baseline experiments P αs is given by P αs = sin 2 2θ αα sin 2 1.27 (23) Here, sin 2 2θ αα = 4|U α4 | 2 (1 − |U α4 | 2 ) (U is the 4 × 4 lepton mixing matrix) and ∆m 2 41 = m 2 4 − m 2 1 is the sterile-active neutrino mass-squared difference. The oscillation probability for active states by P αβ = sin 4 θ αβ sin 2 1.27 where sin 4 θ αβ = 4|U α4 | 2 |U β4 | 2 . The recoil spectrum induced by neutrinos of flavor α is then given by To a fairly good approximation P αβ can be neglected due to the higher-order active-sterile mixing suppression. N T is the number of target nuclei and dσ/dE r is the SM cross section in Eq. (1). The total number of counts in the k th bin is obtained from Eq. (25) according to Experimental information on R can then be mapped into sin 2 θ i j − ∆m 2 41 planes. NGI follows the same approach as NSI, but includes all possible Lorentz-invariant structures [34]. It was introduced in the analysis of neutrino propagation in matter in Ref. [35], studied in the context of CEνNS physics in Ref. [36] and in the light of COHERENT data in Ref. [10]. Dropping flavor indices, the most general Lagrangian reads where Γ a = {I, iγ 5 , γ µ , γ µ γ 5 , σ µν }, with σ µν = i[γ µ , γ ν ]/2. As in the NSI case, some of these couplings lead to spin-suppressed interactions which we do not consider. Relevant couplings therefore include all Lorentz structures for the neutrino bilinear and only scalar, vector and tensor structures for the quark currents. For the NGI analysis, we consider only one Lorentz structure at a time and assume the C and D parameters to be real. We may therefore consider the individual cross sections. Assuming a spin-1/2 nuclear ground state and ne- For the scalar and tensor cases the SM cross section has to be added. In the vector case ξ 2 V includes the SM contribution which interferes with the BSM vector piece. The definition of the different parameters in Eq. (28) can be found in Appendix B.

B. Impact of neutron form factor uncertainties
We calculate the impact of uncertainties in the neutron rms radii for CsI, Ge and Xe in the presence of NSI. To do so, we take as "experimental" input the number of events predicted by the SM assuming r n rms = r p = ∑ i r p i rms X i = 4.06 fm for all germanium isotopes and 4.79 fm for all xenon isotopes. Here r p i rms is the rms radius of the proton distribution of the i th isotope with abundance X i . We proceed as we have done in Section IV, i.e., for CsI we take into account the Cs and I contributions, while for Ge and Xe the contributions for each isotope according to Eq. (18). For all three cases we assume four years of data taking. For our analysis we define the χ 2 function, where α is a nuisance parameter that accounts for uncertainties in the signal rate, N meas i is the number of simulated events in the i th bin, and N BSM i is the number of predicted events in the BSM scenario (which depend on the set of parameters P ). The statistical uncertainty in the simulated data is σ i = N meas i and we have fixed σ α = 0.28 in all cases [1].
Assuming ε q µµ ⊂ [−1, 1] we determine the 90% C.L. exclusion regions in two cases, r n rms = r p and r n rms = 5.29 fm. We find that the CsI and Ge detectors are rather insensitive to the choice of the neutron rms radius; the resulting 90% C.L. regions barely change with r n rms . For Xe, the result is quite different. Changing r n rms has a strong impact on the available parameter space. This can be seen from the left panel of Fig. 5, where the diagonal bands are obtained in the case r n rms = r p , while the purple regions are obtained for r n rms = 5.29 fm. This result is as expected. Firstly, the LXe detector has a larger target mass (about a factor 6.5 larger compared to the CsI and Ge detectors), so for a common data taking time the accumulated statistics in the LXe detector is larger. Secondly, the number of events expected in the NSI scenario with r n rms = r p reproduces the simulated data better than with r n rms = 5.29 fm since the data are simulated with r n rms = r p . In what follows we consider only COHERENT phase-III. Note, however, that increasing the exposure for the CsI or Ge detectors will change the situation. In doing so these detectors will become-as the LXe detector-sensitive to uncertainties in the neutron rms radii.
To determine the extent to which neutrino NSI can be distinguished from the SM signal including its neutron form factor uncertainties, we calculate the number of events assuming ε d µµ = 0 and ε u µµ ⊂ [−0.088, 0.37]. These values correspond to the 90% C.L. range obtained from global fits to neutrino oscillation data including COHERENT (CsI phase) data without accounting for form factor uncertainties [6]. The result is shown in the right panel of Fig. 5. The NSI (purple) histograms are obtained by fixing r n rms = r p . The SM histograms (orange) are instead obtained by fixing r n rms = r p (upper boundary) and r n rms = 5.29 fm (lower boundary) and determines the SM expectation within the form factor uncertainties. Clearly, the SM expectation with form factor uncertainties lies within the NSI expectation for ε u µµ between −0.08 and 0.2 with a fixed form factor. There are various ranges of NSI couplings that will produce signals that cannot be disentangled from the SM signal. This will persist unless uncertainties on the neutron rms radii are reduced. We have chosen ε u µµ to stress this point, although results for ε d µµ , ε q ee , ε q ττ and ε q eτ , will lead to the same conclusion. Needless to say, allowing for multiple nonzero NSI parameters will further complicate the ability to discriminate new physics from the SM.
For sterile neutrinos we calculate the 90% C.L. exclusion regions in the sin 2 2θ 24 − ∆m 2 41 plane. The results are shown in the left panel of Fig. 6. The contours are obtained for N BSM calculated for r n rms = r p (orange contour) and r n rms = 4.9 fm (purple contour). We assumed these values for all xenon isotopes and fix sin 2 θ 14 = 0.01 (best-fit value from a global fit to ν e andν e disappearance data [37]) and θ 34 = 0. This result demonstrates that the available regions in parameter space have a strong dependence on the neutron rms radii. A ∼ 2% change in r n rms is sufficient to significantly modify the results of the parameter fit. For r n rms = 5.29 fm, the fit to the simulated data is extremely poor.
It is clear that a more precise treatment of sterile neutrino effects should include neutron form factor uncertainties, otherwise one might end up misidentifying SM uncertainties with these effects. To show this might be the case, we calculate the number of events for sterile neutrino parameters fixed as in the previous calculation and for ∆m 2 41 = 1.3 eV 2 , r n rms = r p , sin 2 2θ 24 = 0 and sin 2 2θ 24 = 1. We then compare the resulting histograms (in purple) with the SM prediction including uncertainties (in orange); see the right panel of Fig. 6. The overlapping spectra show that an identification of the new effects is not readily possible.
Finally, we turn to the discussion of the impact of the uncertainties on the sensitivity to neutrino NGI. The results are shown in Fig. 7. We have proceeded in the same way as that for the NSI analysis, fixing r n rms = r p to generate the simulated data, and analyzing the results for two cases with r n rms = r p and r n rms = 5.29 fm. For scalar interactions we assume D q P = 0, while for vector interactions, D q A = 0. Results in the case r n rms = r p are quite similar to those found in Ref. [10] and largely depart from them for r n rms = 5.29 fm for reasons similar to that for NSI. N meas has been simulated assuming r n rms = r p and so in the presence of NGI a better fit is found for the first sample. With r n rms = 5.29 fm there is little room for new interactions since the mismatch between the neutron rms radii induces substantial departure from the simulated data. It is clear that scalar NGI interactions are far more sensitive to neutron form factor uncertaintes; changes of ∼ 20% in r n rms almost entirely close the parameter space. Vector and tensor NGI are not as sensitive as scalar NGI, but depending on the value of r n rms large portions of parameter space are allowed or disfavored.

VI. CONCLUSIONS
We have quantified the uncertainties on the SM CEνNS cross section. They are driven by the neutron form factor through its dependence on the rms radius of the neutron density distribution, r n rms . To quantify these uncertainties we assumed that r n rms ranges between the rms radius of the proton charge distribution, r p rms , of the corresponding nuclide and r n rms ( 208 Pb) (which has been precisely measured by PREX). Under this assumption we evaluated the size of the uncertainties for 133 Cs, 127 I, germanium, xenon and argon-choices motivated by COHERENT phases I-III and DarkSide-20kusing three form factor parameterizations: Helm, Fourier transform of the symmeterized Fermi function and Klein-Nystrand. We showed that form factor uncertainties: (i) are relevant for q 20 MeV, and so are negligible if the CEνNS process is induced by either reactor or solar neutrinos, (ii) have percentage uncertainties that reach their maximum at q 65 MeV and are order 8% − 9% for Cs, I, Xe and Ge, and 11% for Ar; the lighter the isotope the larger the percentage uncertainty, (iii) are basically independent of the parametrization used.
We studied the impact of the uncertainties on the SM prediction for COHERENT, diffuse supernova neutrino background, and sub-GeV atmospheric neutrinos. For COHER-ENT, assuming n POT as in [1], we found that the SM prediction is subject to relative uncertainties that are never below 2.5% in germanium and 5% in xenon and that can be as large as 23% or 52% (depending on the recoil energy) in the germanium and xenon detector, respectively. For DSNB neutrinos we find that the relative uncertainties range from 4% to 7%, while for sub-GeV neutrinos the uncertainties are up to 16%. These results demonstrate that in the absence of precise measurements of r n rms , SM predictions of the CEνNS rate involve uncertainties that challenge the interpretation of data. This is especially true for future measurements with small experimental systematic uncertainties.
We also quantified the impact of the neutron form factor uncertainties on the sensitivity to new physics. We considered three scenarios: neutrino NSI, sterile neutrinos in the 3+1 scheme, and NGI. We showed that the variation of r n rms has a profound effect on these new physics searches.
Finally, it is worth pointing out that the uncertainties we have derived here also apply to DM direct detection searches, provided the WIMP-nucleus interactions are spin independent. In WIMP scenarios with vector, scalar and tensor mediators, the direct detection rate will involve uncertainties comparable to those we have derived.
In core-collapse supernova, neutrinos of all flavors are emitted and each flavor carries about the same fraction of the total energy, E ν 3 × 10 53 erg. Their spectra are approximately thermal with temperatures obeying Tν e < T ν e < T ν x (ν x = ν µ , ν τ ,ν µ ,ν τ ). We have taken a Fermi-Dirac distribution with zero chemical potential for all flavors,