Determining the nuclear neutron distribution from Coherent Elastic neutrino-Nucleus Scattering: current results and future prospects

Coherent elastic neutrino-nucleus scattering (CE$\nu$NS), a process recently measured for the first time at ORNL's Spallation Neutron Source, is directly sensitive to the weak form factor of the nucleus. The European Spallation Source (ESS), presently under construction, will generate the most intense pulsed neutrino flux suitable for the detection of CE$\nu$NS. In this paper we quantify its potential to determine the root mean square radius of the point-neutron distribution, for a variety of target nuclei and a suite of detectors. To put our results in context we also derive, for the first time, a constraint on this parameter from the analysis of the energy and timing data of the CsI detector at the COHERENT experiment.

Since CEνNS is sensitive to the weak form factor, dominated by the coupling to neutrons, this process can probe the distribution of neutrons in nuclei. This is precious information to complement proton densities accessible with elastic electron scattering [42,43]. At present the most direct measurement of a neutron distribution comes from parity-violating electron scattering in 208 Pb [44], also sensitive to the weak form factor. Alternative measurements rely on nuclear [45][46][47][48][49] or electromagnetic [50] reactions which probe both neutron and proton distributions, but they lean on model-dependent analyses (with uncertainties that are difficult to quantify). Likewise, atomic parity violation experiments are also sensitive to the nuclear neutron distribution, but they are subject to model-dependent uncertainties from atomic many-body calculations [51][52][53][54].
Therefore, CEνNS can shed light on neutron distributions in nuclei, in particular their neutron radii. The difference between the radii of neutron and proton distributions is called neutron skin thickness, or just neutron skin. Its understanding impacts the limits of existence [55] and size [56] of atomic nuclei, and serves as an important test of firstprinciples nuclear calculations [49,57]. Beyond the structure of nuclei, the neutron skin can be related to the energy needed to form isospin asymmetric nuclear matter-the symmetry energy-and its variation with the nuclear density [58][59][60][61]. These are key properties of the equation of state of neutron-rich matter, which determines the size and structure of neutron stars [62,63].
The opportunity to complement and improve over the first measurements of CEνNS using the upcoming European Spallation Source (ESS) in Lund (Sweden) has been recently highlighted in Ref. [64]. The ESS will generate the most intense neutron beams for multidisciplinary science, and an order of magnitude increase in neutrino flux with respect to the SNS. Using novel detector technologies stemming from recent advances in dark matter and neutrinoless double beta-decay experiments, the proposed CEνNS@ESS would be able to maximally profit from the much higher statistics available at the ESS.
Reference [64] presented the physics potential to constrain new physics scenarios in the neutrino sector using CEνNS@ESS. In the present work, on the other hand, the physics scope is very different: we focus the study of nuclear structure with the aim to determine the root mean square (rms) radius of the neutron distribution (R pt n ) in the target nucleus, for several detector materials (CsI, Xe, and Ge). In doing so the treatment of systematic uncertainties is particularly relevant. In this respect, technically, an important difference with respect to Ref. [64] is that, in the present work, we also consider the effect of the energy scale uncertainty, which, as we will show, has a non-negligible impact on the determination of R pt n .
In order to put our results into a larger context, we also derive the present bounds on R pt n for CsI obtained from the analysis of the COHERENT CsI data using both energy and timing information. While other authors have studied the bounds on R pt n using energy infor-mation alone, to our knowledge this is the first time that energy and timing information are used to derive a bound on R pt n . As we will see, this brings additional synergies onto the table and helps to improve the constraint with respect to the case where only energy information is used when the analysis is performed under the same assumptions. Furthermore, in our fit, we also implement an improved background model and several choices of the quenching factor (that gives the relationship between the number of photoelectrons detected and the nuclear recoil energy of a given event). In particular, following Ref. [17] we make use of the new parametrization of the TUNL data which results on a new quenching factor function never considered in the extraction of R pt n . We will show how the inclusion of these effects leads to a spread in the extracted value of R pt n beyond those previously considered in the literature. Futhermore our results illustrate that current theoretical determinations of the neutron radius may favour some values of the quenching factor.
The article is structured as follows. In Sec. 2 we introduce our notation and describe the methodology used in our numerical fits and simulations. Our results are presented in Sec. 3, both for current bounds obtained using the COHERENT CsI data (Sec. 3.2) and for future sensitivities expected at the ESS (Sec. 3.3). Finally, we summarize and conclude in Sec. 4.

Notation and framework
The differential cross section for CEνNS , for a neutrino with incident energy E ν on a nucleus of mass M , can be generically written as [65] where G F is the Fermi constant, T is the recoil energy of the nucleus (kinematically restricted to the interval T ∈ 0, 2E 2 ν /(M + 2E ν ) ) and F V,A are the vector and axial form factors of the nucleus which are functions of the squared tree-momentum transfer Q 2 ≡ | q| 2 = 2M T + O T M . G V and G A are the effective vector and axial couplings for the effective weak current of the neutrino-nucleon interaction. The axial part of the cross section is sensitive to the distribution of nucleon spins in the nucleus. Because of the attractive nuclear pairing interaction, the sum over nucleon spins cancels to a large extent [66], so that axial terms are typically smaller by a factor 1/N 2 . Furthermore, nuclei with even number of protons and neutrons have zero spin, so that axial terms vanish. With this, the differential cross section for CEνNS can be conveniently written as where we have introduced the weak charge form factor of the nucleus F W ≡ F V , as is commonly done in the literature. For a nucleus with Z protons and N neutrons, G V can SM (tree level) SM expected value Table 1. SM values of the up-(u) and down-quark (d) couplings to the Z boson, for left-and right-handed particles (L and R, respectively); as well as the vector couplings for protons (p) and neutrons (n). Here, s 2 W ≡ sin 2 θ W , where θ W is the weak mixing angle. The expressions on the left column correspond to the tree-level couplings in the SM. The values on the right column are taken from Ref. [67] and include propagator as well as vertex and box corrections, as detailed in Ref. [68].
be written as a linear combination of the fundamental NC couplings of the quarks and neutrinos: In these expressions, the left-handed neutrino NC couplings have already been substituted into G V , while the left-and right-handed quark NC couplings f qL and f qR (q = u, d) as well as the vector proton and neutron NC couplings g p V and g n V are summarized in Tab. 1 for convenience.
Neglecting relativistic Darwin-Foldy [69] and spin-orbit [70] corrections, typically below 0.1% for the relevant Q 2 values in CEνNS [70][71][72][73], the weak form factor of the nucleus can be written as where r 2 p and r 2 n are the squared charge radii for the proton and neutron, respectively. The terms in parentheses in Eq. (2.4) represent the lowest-order nucleon form factors, while F M p,n (Q 2 ) stand for the nuclear structure factors. In particular, F M p (Q 2 ) and F M n (Q 2 ) are the spin-independent proton and neutron structure factors, respectively (see, e.g., Refs. [72,73] for details). These encode the response of the nucleus to the interaction, taking into account that the scattering takes place with a complex many-body system. Their normalization is such that F M p (0) = Z and F M n (0) = N .
Several phenomenological parametrizations exist in the literature for the nuclear structure factors, such as the Helm [74], symmetrized Fermi [75] and the Klein-Nystrand [76] parametrizations, among others. In our simulations, for concreteness, we use the Helm parametrization [74] where j 1 is the first order spherical Bessel function, R 2 0 ≡ 5 3 R 2 W − 5s 2 and s is set at s = 0.9 fm [77]. Nevertheless, given the low momentum transfers involved in CEνNS, it is enough to characterize the structure factors using the first moment of the distribution in Q 2 . For the nuclear structure factors F M i (where i = p, n) this corresponds to the so-called point-proton and point-neutron distribution radii 1 [40,70,71,78] Likewise the weak form factor can be expressed in terms of the weak radius R W : where, according to our normalization of the form factors, F W (0) = 1. The weak radius can be then expressed in terms of R pt p , R pt n , r 2 p and r 2 n , as Equation (2.8) can be simplified by introducing the charge radius R 2 ch , which is precisely determined from elastic electron scattering. Defined in terms of the charge form factor as in Eq. (2.6), with same approximations as Eq. (2.8) the charge radius reads [79] This yields (2.10) From Eq. (2.10) we directly read that, in addition to the information accessible in electromagnetic scattering experiments, CEνNS provides independent information on the difference between the rms radii of the neutron and the proton distributions, the neutron skin. Knowledge of the neutron skin is important to learn about nuclear structure and test nuclear models [32,49,[55][56][57]60]. In addition, it can constrain the equation of state of neutron-rich matter [58][59][60][61], a key ingredient for the structure of neutron stars [62,63].

CsI
Xe Ge  Table 2. Charge radius (R ch ) and abundances (%) used for the different isotopes considered in this work. While for the CsI detectors we assume a 50% number abundance of Cs and 50% of I, for Xe and Ge we use their natural isotope abundances from Ref. [82]. The values for the charge radii, taken from Ref. [43], are measured to a few per mil precision, and therefore we ignore their error bars in our calculations.
From Eqs. (2.7) and (2.10) it is easy to see that larger (smaller) values of R pt n tend to suppress (enhance) the number of events with large momentum transfer. Since Q 2 2M T this, in turn, affects the shape of the event distribution as a function of the recoil energy of the nucleus, leading to a similar suppression/enhancement of the tail.
In our simulations, we fit the data (either real data from COHERENT, or simulated data for the ESS) to extract the weak radius, and use Eq. (2.10) to obtain the rms radius of the point-neutron distribution R pt n . In doing so, we use as inputs the tabulated nuclear charge radii R ch from Ref. [43], together with the Particle Data Group (PDG) values for the proton and neutron charge radii 2 : r 2 p = 0.70706 fm 2 and r 2 n −0.1161 fm 2 [67]. For convenience, the values of R ch used in our calculations are summarized in Tab. 2. Finally note that, in the case of CsI, we assume the same rms radius for the neutron distributions of Cs and I, because only one combined number can be extracted from CEνNS . This assumption seems reasonable given the sensitivity of current experiments and the similar values of Z and N for Cs and I.

Signal and background event rates for a CEνNS experiment
At both the SNS and the ESS, the neutrino flux is predominantly produced from pion decay at rest: while π − get rapidly absorbed by the nuclei after being produced, the π + eventually decay at rest into π + → µ + ν µ . The prompt (monochromatic) ν µ flux component is followed by a delayed contribution fromν µ and ν e , produced in the decay of the muon µ + →ν µ ν e . The yearly energy spectrum of neutrinos reaching a detector at a distance from the source, summed over all neutrino flavours, reads (3.1) where E ν is the neutrino energy, m µ and m π are the muon and charged pion masses, respectively, N POT is the number of protons on target (PoT) delivered per year, f ν/p is the neutrino yield per PoT (directly related to the pion yield per PoT), δ is the Dirac delta function and θ is the Heaviside function.
At a CEνNS experiment, three main sources of backgrounds should be considered: (i) steady-state (SS) backgrounds (dominated by cosmic ray interactions or by their byproducts inside or in the surroundings of a radio-clean detector); (ii) beam-related backgrounds, produced by neutrons escaping the target and reaching the detector; and (iii) neutrino-induced neutrons, that is, neutrons produced in neutrino interactions inside (or in the surroundings of) the detector. At COHERENT, the last background contribution was determined to be very small [2], and therefore will be also neglected here for the ESS. Beam-related background are also expected to be sufficiently suppressed, see Ref. [3]. While we include their expected contribution in our fits to COHERENT, they will be neglected in our ESS simulations for simplicity. The most relevant background will therefore be the SS contribution. In the case of COHERENT, this is estimated and modeled using data taken when the proton pulse is turned off, as we will discuss in more detail below. At the ESS, the possible contribution from this background can only be estimated from Monte Carlo simulations, which however are not yet available. Therefore, we conservatively assume it to be uniformly distributed in T , reading its normalization from Tab. 1 in Ref. [64].

Present bounds from COHERENT energy and timing data
For the sake of completeness and comparison with the existing literature we start by performing an analysis of the present results from COHERENT experiment using the detailed timing and energy information of their data for CsI. In particular, we analyze the data provided in Ref. [83] for 308.1 live-days of neutrino production, which corresponds to 7.48 GWhr (or ≈ 1.76×10 23 protons on target). The events observed are binned in a two-dimensional grid, using the number of photoelectrons observed (equivalent to the nuclear recoil energy) and the time with respect to the start of the beam pulse. Although COHERENT data was used in Refs. [15,32] to extract the rms neutron radius, to our knowledge this is the first time the timing information is included in the determination of the neutron distribution radius in CsI. As we show below, this has an impact on the results.
We refer to the reader to Ref. [17] for details on our analysis of the COHERENT CsI data, while here we just summarize some of the most relevant details of the fit. In brief, we perform a fit to the data following the COHERENT data release albeit with some modifications. Namely, we study the dependence of the results with respect to the choice of quenching factor (QF) as well as with the implementation of the steady-state background component. In particular we perform the analysis for 3 choices of the QF and its associated uncertainties: (i) the energy-independent QF used in the COHERENT  Figure 1. ∆χ 2 as a function of the rms radius of the point-neutron distribution R pt n for Cs or I (assumed to be equal), for a variety of fits to COHERENT data as indicated by the labels. In all cases shown in the left panel, the QF and background assumptions are those employed in the data release [3]. On the right panel we show the dependence of the results on the assumptions for SS background modeling and QF implementation. For convenience, the vertical line indicates the value of the rms radius for the proton distribution. This is taken as the average between their values for Cs (R pt,Cs data release; (ii) the energy-dependent QF with reduced uncertainties obtained in Ref. [84] (hereafter referred to as "Chicago QF"); and (iii) a new QF parametrization from our own fit to the calibration measurements performed by the TUNL group [2,3] (hereafter referred as "Duke QF"). In what respects to the background treatment, we use two different parametrizations for the temporal behaviour of the SS background: (a) we use the ad-hoc exponential parametrization prescribed in the experimental data release [3] (which however leads to an unexplained mild excess in the first two time bins); and (b) we use our own parametrization of the temporal behaviour of the background, based on a temporal fit to the number of events detected when the neutrino beam is turned off.
Our results for the fit are shown in Fig. 1, where we plot the dependence of the χ 2 on R pt n for the different models for the QF and background implementation used to fit the data. In the left panel we illustrate the dependence of our results on R pt n , for a fit which uses the same QF and SS background parametrization as in the official data release. The different lines show the results obtained using only the total rate, energy information, and/or timing information, as indicated by the labels. Two salient features are identified right away from this panel. First, comparing the dot-dashed and dotted lines, we find that the inclusion of the energy dependence of the data does not lead to a substantial improvement for the determination of R pt n . This is so because of the large systematic uncertainties assumed for the QF employed in the data release, which do not allow to observe a significant variation in the shape of the event distributions as R pt n is varied. In contrast, once timing information is included the χ 2 increases significantly for small values of R pt n , shifting the best-fit towards slightly larger values. This is mainly driven by the small excess observed for the event rates in the first two time bins (with respect to the SM prediction), which is in mild tension with the prediction if we use an exponential fit to model the SS background (as prescribed in the data release, see Ref. [17] for details). Therefore, the excess could in principle be accommodated with a corresponding increase of the prediction for the prompt neutrino contribution. Note that the prompt neutrino flux is characterized by lower energies and therefore leads to lower values of T , while the energies of the delayed component of the flux are larger and contribute more to the high-energy tail of the observed spectrum. Although the value of R pt n affects all neutrino species (and therefore the prompt and delayed signals) in the same way, a large value of R pt n would suppress the tail of the distribution for large nuclear recoils; this, combined with an overall increase of the total normalization (thanks to the large systematic uncertainties affecting the fit) would effectively induce a change in shape for the event distribution, mimicking an enhanced prompt neutrino contribution.
On the other hand, the right panel of Fig. 1 quantifies the effect on the fit due to changes in the QF and the background treatment. First, as forecasted in the above discussion, once the ad-hoc exponential parametrization of the temporal behaviour of the SS background is substituted by a data-based modeling, the prediction is able to fit better the observed data, and the induced excess at low time disappears. As a result the best-fit is slightly shifted to lower values, even when using the same QF implementation and associated systematic error as in the data release. From the comparison between the solid red (Data release t + E) and solid orange (+ new Background) lines one can also see that the size of the confidence regions is also slightly larger for the latter, as it is often the case when the tension between data and background model is eliminated. Using our own fit to the TUNL data (Our fit t + E D) results into a reduction of the uncertainty (as expected, since the size of the systematic error associated to the QF is much smaller in this case) but maintains the best-fit at the same value. Instead, if the Chicago QF is used (Our fit t + E C), the allowed range of R pt n shifts to significantly lower values. The reason is that this QF leads to a reduced number of predicted events (see Fig. 2 in Ref. [17]), which fall somewhat short to explain the observed data. In the fit, this can be partially compensated by allowing for a smaller R pt n , which leads to an enhancement of the number of events.

Future prospects at the ESS
The ESS will soon generate the largest pulsed neutrino flux suitable for the detection CEνNS. In a recent article [64] the potential for particle physics phenomenology was quantified for the ESS, using a series of innovative detector technologies specifically designed to detect very low nuclear recoils, as those expected for CEνNS . Here we summarize the main differences between the ESS and SNS neutrino sources, and describe how we adapt the simulations performed in Ref. [64] to the study of the neutron distribution radius, while we refer the reader to Ref. [64] for additional details.
In what respects the assumed neutrino flux, the ESS is scheduled to reach its design power of 5 MW by 2023. This is to be compared to the nominal 1 MW power of the SNS. Considering the higher proton energies at the ESS (2 GeV, compared to 1 GeV at the SNS), this represents an increase in average proton current at the ESS by a factor of 2.5, for a total of N POT = 2.8 × 10 23 protons on target per calendar year and approximately 5,000 hours of beam delivery per year. An additional increase in the total neutrino flux arises from the larger pion yield per proton at the ESS, due to the much higher proton energies envisioned. Based on a set of simulations performed in Ref. [64] we tentatively adopt a yield of f ν/p = 0.3 neutrinos of each flavor (ν µ ,ν µ , ν e ) per proton for a ESS operating at 2 GeV, which is about 3.2 times larger than the corresponding yield at the SNS. Conversely, the time length of the proton pulse at the ESS is much broader than that of the SNS. This prevents the use of timing information for flavour discrimination. Although this affects significantly the sensitivity to some BSM scenarios, since the value of R pt n affects the cross section for all neutrino flavors in the same way, this will not be so relevant for our physics case at hand. On the other hand, the much longer proton pulse yields a duty factor of 4 × 10 −2 , almost two orders of magnitude larger than the SNS duty factor of 6 × 10 −4 . Therefore, during the beam spill time window a larger number of SS background events will enter the detector. Although these backgrounds can be well-characterized using beam OFF data, their associated statistical uncertainties have to be properly accounted for in the simulations.
We compute the expected sensitivities for three different type of detectors which use three different nuclear targets: (a) a cryogenic undoped CsI scintillator array; (b) a highpressure gaseous Xe chamber; and (c) a low-threshold, multi-kg p-type point contact Ge detector. The rationale behind this choice is that the first CsI detector allows for a direct comparison with the present bound from COHERENT presented in the previous section, since both use the same target material; while the other two detectors use a heavy (Xe) and medium-mass (Ge) nuclear targets, allowing us to quantify the impact on the sensitivity due to the choice of target nucleus. Finally, since the optimal detector location for the ESS has not been identified yet, in all our simulations we assume the detector distance to be = 20 m as a reasonable benchmark value. In all cases we assume an exposure of 3 years, corresponding to 8.4 × 10 23 total PoT. Table 3 summarizes the main characteristics relevant to the simulations, while we refer the interested reader to Ref. [64] for detailed descriptions and characterizations of these detectors.
Our sensitivity calculations for the ESS use event distributions that are binned in nuclear recoil energy. At a CEνNS experiment, the typical observable is usually the number of detected photoelectrons (PE) in an event. The number of PE is related to the nuclear recoil by the QF, which is determined from experimental calibration measurements. For COHERENT, since the data are provided in terms of the number of PE we have performed the analysis using that variable. Conversely, since for the ESS we are dealing with a future proposal, we have decided to bin the data in recoil energy instead. It should be stressed out that, once the QF is known, both variables are completely equivalent and the sensitivity analysis would give the same results regardless of the variable used to bin the data. Also note that, while we do not use a specific QF in our simulations for the ESS, we do consider an associated systematic error, as described below.
Within each bin i with reconstructed nuclear recoil energy T rec ∈ [T i , T i+1 ], the expected  Table 3. Summary of detector properties and background rates used in our sensitivity calculations. From left to right, the different columns indicate the target nucleus, total detector mass, steadystate (SS) background rates, nuclear recoil detection threshold in keV, energy resolution at threshold, maximum recoil energy considered, and the total number of bins used in the simulation (see text for details). Backgrounds rates are listed in counts per keV, kg and day (ckkd), before applying the 4 × 10 −2 reduction due to the ESS duty factor. We conservatively adopt a flat SS background estimated at T th (where it is typically largest) for the whole energy range.
number of events per year N i can be computed as is the number of target nuclei j in the detector. While for the CsI detector we assume a 50% number abundance of Cs and 50% of I, for Xe and Ge we use their natural isotope abundances as taken from Ref. [82], see Tab. 2.
In Eq. (3.2), R(T rec , T ) is the the energy reconstruction function. We assume it to be a Gaussian with a width that depends on the recoil energy as: σ(T ) = σ 0 T /T th , where σ 0 is the energy resolution at the detection threshold (T th ), see Tab. 3. For each detector, the recoil energy bin sizes are chosen so that the width of each bin is twice the energy resolution at its center. We consider all the kinematically available range (determined by the condition T 2E 2 ν /M ), for all detector configurations. For convenience, the total number of bins is provided in the last column of the Tab. 3.
In order to determine the sensitivity to R pt n we first simulate the future data as the expected number of events per bin,N i , for a given detector and for an assumed true value for the rms neutron radius, R pt,true n . For sake of concreteness, we assume R pt,true n = 1.05R pt p ; however, we have numerically checked that our results do not depend significantly on the chosen value for R pt,true n as long as it lies within the range ∼ ±25% R pt p . Once the event distribution for the data (N i ) has been simulated, and the expected event rates (N i ) have been computed as a function of R pt n , a binned χ 2 is built. Systematic uncertainties are implemented using the pull method. In doing this, a set of nuisance parameters η j is included in the fit, each of them with an associated prior σ j . A penalty term is then added to the χ 2 for each source of systematic errors, and the minimum of the χ 2 is determined after marginalization over all the nuisance parameters included in the fit. We consider three different types of systematic uncertainties: 1. An error on the total signal normalization. To implement it, we substitute N POT → N POT (1 + η norm σ norm ) in Eq. (3.2). We set σ norm = 0.1 [64].
2. An error on the total background normalization. To implement it, we substitute N bkg . We set σ bkg = 0.05 [64].
3. A systematic error affecting the energy scale (ES), directly related to the QF uncertainty. To implement it, we substitute T rec → T rec (1 + η ES σ ES (T rec )) in Eq. (3.2).
Here, σ ES is the prior ES uncertainty, which may depend on the reconstructed recoil energy. Since we find that this uncertainty is the one that has the largest impact on our results (see below for details), we will consider three different assumptions for σ ES : negligible; an energy-independent 5% uncertainty; and an energy-dependent uncertainty which decreases linearly from 5% at T rec,th to 1% at T rec,max .
With all these ingredients, our χ 2 function reads (3.3) whereN i are the assumed data event rates at face value, i.e.,N i = N i (R pt n = R pt,true n , η j = 0). Figure 2 illustrates the impact on the sensitivity due to the ES uncertainty. In both panels, the black dots show the simulated data for the CsI detector, for our assumed value of R pt,true n . The colored histograms, on the other hand, show the predicted spectra for different test values of R pt n (R pt n < R pt,true n in the left panel, while R pt n > R pt,true n in the right panel), and for different assumptions on the ES systematic error, after setting the nuisance parameters at the values which give the best possible fit to the simulated data points. The dotted green histograms have been obtained without an ES uncertainty, and therefore show directly the impact of R pt n on the event distributions: as expected from the analytic expressions in Sec. 2, an excess of events is observed for R pt n ≤ R pt,true n (left panel), particularly relevant at high recoils, while the right panel shows the opposite behaviour. Once an ES uncertainty is added, the fit will try to vary the nuisance parameters in the χ 2 to find a better fit to the data. As can be seen, the inclusion of this pull term allows to induce apparent changes in the observed spectrum, which mimic the impact of a different value of R pt n in the fit. This is better appreciated in the lower panels, which show the relative difference in the event rates per bin (with error bars showing the corresponding uncertainties). As shown in the lower panels, the inclusion of an energy-independent ES uncertainty significantly relaxes the tension in the fit in both cases.
Finally Fig. 3 shows the resulting χ 2 (R pt n ) for our ESS simulations, as a function of R pt n /R pt,true n for the three different detectors as well as for the different assumptions of the systematic uncertainties considered. By construction, the χ 2 (R pt n ) has its minimum at zero for R pt n /R pt,true the dotted lines in each panel and the dashed/solid lines we immediately observe that the inclusion of an ES uncertainty, even if smaller than the overall normalization uncertainty, significantly spoils the sensitivity. This is expected, since the main effect of R pt n is the distortion of the tail of the recoil energy spectrum as described in Sec. 2 (see also Fig. 2), an effect that can be mimicked by an ES uncertainty.
From the comparison between the different panels in Fig. 3 we can also see how the sensitivity is substantially worsened for lighter nuclei. This is expected a priori since the signal statistics for CEνNS grows quadratically with the number of neutrons in the target nucleus, see Eqs. (2.2) and (2.8). Moreover, note that the values of R pt p,n (or, equivalently, R W ) characterize the nuclear size, and therefore are significantly smaller for Ge than for Xe or CsI. Thus, from Eq. (2.7) one can see that for smaller nuclei the dependence of the form factor with Q 2 will be very mild in these cases. Since the sensitivity to R pt n comes precisely from the observation of a change in the shape of the energy distribution, this automatically translates into a worse sensitivity for smaller nuclei. In other words: the sensitivity to R pt n comes from the fact that the wavelength of the neutrino is of the order of the size of the nucleus, so it is sensitive to the nuclear structure and, in particular, to the size of the neutron distribution. However, as the size of the nucleus decreases larger neutrino energies are required to probe the size of the nucleon distribution. This poses a challenge for CEνNS, for which low momentum-transfers are required to maintain coherence.  Figure 3. ∆χ 2 as a function of the rms radius of the neutron distribution R pt n for three detectors proposed for a CEνNS experiment at the ESS: a cryogenic undoped CsI scintillator array (left panel) a high-pressure gaseous xenon chamber (central panel), and a low-threshold, multi-kg p-type point contact germanium detector (right panel). In all cases an exposure of 3 years is assumed, together with a 10% signal normalization uncertainty and a 5% background normalization uncertainty. The three curves in each panel show our results for different assumptions for the systematic uncertainty in the energy scale reconstruction: negligible (dotted lines), energy-dependent, varying from 5% at T rec,th to 1% at T rec,max (dashed lines), and 5% energy-independent (solid lines).

Summary and Conclusions
Coherent elastic neutrino-nucleus scattering (CEνNS ) probes the weak form factor of the nucleus and provides us with information of its weak charge distribution. In particular, from CEνNS data it is possible to determine the rms radius of the neutron distribution of the target nuclei. The difference with the rms proton distribution-the neutron skin thickness-is crucial in a broad spectrum of nuclear physics and astrophysics, including nuclear structure, nuclear matter, and neutron stars.
First, as reference, and for comparison with the existing literature we have performed an analysis of the present results from COHERENT experiment using the detailed timing and energy information of their data for CsI. We found that COHERENT provides a 1σ uncertainty on the rms radius of the neutron distribution which ranges from 11% to 16%, with an additional 11% uncertainty that depends on the choice of QF and SS background modeling. In particular, we get the following determination of R pt n (at 1σ) from the different analyses: This is to be compared with the corresponding values for the proton distributions for Figure 4. Compilation of our results for the present determination and future sensitivity of the neutron skin thickness from CEνNS experiments for the isospin asymmetric nucleus considered in this work. For clarity, the red and blue points corresponding to two analyses of COHERENT have been slightly displaced horizontally. For comparison we show the result for 208 Pb derived from the PREX measurement with the analysis in Ref. [71], and the best-fit to the indirect determination from antiprotonic atom x-ray data for a variety of nuclei (dashed line), taken from Ref. [48]. The hollow markers represent the range of predictions in a variety of models, see text for details. For concreteness, the expected results from the ESS have been centered on the dashed line, and correspond to our analysis with an energy-dependent ES uncertainty ranging from 1% to 5%.
Cs, R pt,Cs p = 4.75 fm, and I, R pt,I p = 4.70 fm. In Fig. 4 we plot the corresponding results as the inferred neutron skin thickness by subtracting the average R pt,CsI p = 4.725 fm. Despite the relatively poor precision, it is important to stress that these are the only direct probes available for this observable. Even within present uncertainties, it is interesting to point out that the results from the data release (original background) and our fit with the Duke QF lead to values of the neutron skin that exceed those predicted by theoretical nuclear models [32,85]. In contrast, the analysis with the Chicago QF yields a thinner neutron skin, consistent with calculations. However, note that the Chicago QF 1σ errors reach negative values for the inferred neutron skin thickness, an unexpected result for a neutronrich nucleus [86] not predicted by any nuclear model [32,55,57,60,87].
For comparison, Fig. 4 shows the neutron skin thickness for 208 Pb, derived from the analysis in Ref. [71] of the results of the Lead Radius Experiment (PREX) experiment [44] on parity-violation in electron scattering. Let us stress that although parity-violation is a weak measurement (and therefore it does probe the neutron distribution of the nucleus), PREX does so only at a fixed momentum transfer. The extraction of the point-neutron radius and the corresponding skin thickness from the asymmetry cross-section measured by PREX contains certain model dependence, as can be observed by the slightly different results presented by Refs. [44] and [71]. In this figure we also show the best-fit for the neutron skin obtained from its indirect determination from antiprotonic X-ray data for a variety of nuclei (dashed line), taken from Ref. [48]. Finally, the hollow markers show the predictions for a variety of nuclear models as extracted from Refs. [32,73,85,87].
Regarding the future sensitivity at a CEνNS experiment using the ESS as the neutrino source we have focused on the proposed detectors using three different detector materials: CsI, Xe and Ge. With the large statistics expected, variations of R pt n of O(few %) can lead to observable distortions in the recoil energy spectrum provided it can be measured with enough precision. We have explored this effect and obtained the following expected precision for the determination of R pt n at 1σ, for the different detectors and three choices of energy scale (ES) systematic uncertainties, namely, negligible ES uncertainty (5% to 1% ES uncertainty) [5% ES uncertainty]: Our results for the expected precision at the ESS are also shown in Fig. 4, for the analysis with (5% to 1% ES uncertainty). While the uncertainties in Eq. (4.5) lead to a range of values of the neutron skin thickness which exceeds the typical spread of nuclear structure results using different models, ∼ (0.1 − 0.3) fm [32,60,87], Fig. 4 shows that CEνNS can provide useful constraints on the neutron skin of heavy, neutron-rich nuclei. In fact CEνNS gives the most direct measurement of the neutron skin, with nuclear-model independent error bars, bringing the opportunity to measure it for different nuclei relatively easily. Combined with future improvements beyond the ESS, these advantages place CEνNS in a position towards determining unambiguously the size of the neutron distribution with model-independent uncertainties.