A global analysis strategy to resolve neutrino NSI degeneracies with scattering and oscillation data

Neutrino non-standard interactions (NSI) with the first generation of standard model fermions can span a parameter space of large dimension and exhibit degeneracies that cannot be broken by a single class of experiment. Oscillation experiments, together with neutrino scattering experiments, can merge their observations into a highly informational dataset to combat this problem. We consider combining neutrino-electron and neutrino-nucleus scattering data from the Borexino and COHERENT experiments, including a projection for the upcoming coherent neutrino scattering measurement at the CENNS-10 liquid argon detector. We extend the reach of these datasets over the NSI parameter space with projections for neutrino scattering at a future multi-ton scale dark matter detector and future oscillation measurements from atmospheric neutrinos at the Deep Underground Neutrino Experiment (DUNE). In order to perform this global analysis, we adopt a novel approach using the copula method, utilized to combine posterior information from different experiments with a large, generalized set of NSI parameters. We find that the contributions from DUNE and a dark matter detector to the Borexino and COHERENT fits can improve constraints on the electron and quark NSI parameters by up to a factor of 2 to 3, even when relatively many NSI parameters are left free to vary in the analysis.


Introduction
Non-standard neutrino interactions (NSI) are a popular effective field theory framework for exploring new physics beyond the standard model (BSM) in the neutrino sector [1][2][3]. In the context of neutrino scattering experiments and neutrino oscillations, in the limit where any new gauge fields that mediate NSI are much heavier than the characteristic momentum transfer q 2 , they are a convenient expression of the effective operators that arise in BSM extensions. NSI have been considered in many contexts and for a variety of neutral current and charged current operators. We limit the scope of this study to dimension-6 neutral-current (NC) vector NSI among the first-generation of SM fermions. They are described by the effective Lagrangian Here the fermion indices are f = e, u, d and flavor indices α, β = e, µ, τ . P L and P R are the left and right projection operators, respectively. The effective dimension-6 operators in Eq. 1.1 arise from some fundamental renormalizable theory [4] where the NSI parameters f,L , f,R are taken as proxies for the new propagators multiplied by couplings in the M 2 >> q 2 limit. Experiments that are sensitive to different interaction channels and to different neutrino energies constrain fermion, flavor indicies, and projection operators. For example, solar neutrino experiments sensitive to neutrino-electron scattering are primarily sensitive to electron-type NSI, and also place the most stringent bounds on right-handed NSI.
In previous works that performed statistical analyses of NSI, it has been common practice to either consider a large family of NSI but only vary one or two of them at a time in the likelihood fit (Ref. [5], for example), or reparameterize the NSI down to a more phenomenological and pragmatically manageable subset based on model assumptions (for example, in Refs. [6][7][8]). This is usually done for (i) the sake of model simplicity and (ii) computational limitations with regard to the dimensionality of the fit. However, in this work we are motivated to instead take an approach which is substantially more model-independent and generalized to more degrees of freedom.
Regarding (i) we note that a scenario in which more than two NSI are nonzero at once, albeit complex, have no good reason to be prohibited by nature. For those readers that may be interested in how NSI studies can guide model-building in the neutrino sector, a larger NSI parameter space is warranted to provide generalized constraints. Additionally, degeneracies among the NSI parameters arise due to transformations that leave the oscillation Hamiltonian and scattering cross sections invariant. The full space of these degeneracies as they show up in a likelihood analysis are not fully explored if only a small subset of NSI parameters are activated. Therefore, to explore this large-dimension scenario, we aim to perform a global analysis with all real-valued NSI in Eq. 1.1 nonzero. A model-independent analysis of NSI of this breadth has not been performed to date.
To address issue (ii), we have developed a new statistical technique by virtue of divide-and-conquer which allows one to perform a large-dimensional analysis with a variety of experimental data that are sensitive to different linear combinations of the NSI parameters. The tool in question which allows us to pursue this study without technological barriers is the copula, a statistical object popularized in other data-driven fields but which is quite novel to particle physics. We will discuss this technique in detail in Section 3.
While working in the context of a many-parameter NSI study, the degeneracies that present themselves in physical observables motivate a specific combination of experimental data that can break such degeneracies. While looking forward to the plentiful source of neutrino oscillation data at DUNE, we also raise awareness that DUNE's excellent projected sensitivity will only be indirectly sensitive to the electron, u and d quark NSI via their linear combination that enters into the matter potential of the oscillation Hamiltonian. We therefore recognize the need to augment DUNE's future oscillation measurements with neutrino scattering data, namely those from the COHERENT and Borexino experiments, which have more direct access to these NSI. Additionally, as the next generation of multi-ton scale dark matter detectors will be sensitive to neutrino interactions, they also provide a means to study NSI. In particular, natural neutrino sources such as the solar and atmospheric neutrino fluxes contain τ -flavor neutrinos which can complement the ν τ -deficient neutrino source at COHERENT. In this work we envision a unified experimental dataset comprised of neutrino oscillation and scattering data at COHERENT and Borexino, joined with future projections for DUNE and a ton-scale dark matter detector, to carry out a generic, multi-dimensional NSI analysis.
The paper is organized as follows. In Section 2 we break down the varieties of degeneracies among the NSI parameters and discuss how these degeneracies may complement each other in a global analysis. In Section 3 we demonstrate our global analysis strategy and in Section 4 we briefly outline our analysis methods for each experiment under consideration. Finally in Section 5 we present and discuss the posterior distributions of all the NSI parameters included in the analysis and in Section 6 we conclude.

Degeneracies
The phenomenology of neutrinos scattering with SM fermions in the first generation exhibit several experimental degeneracies, or transformations in the NSI parameter space that leave a physical observable such as a cross-section or a Hamiltonian invariant. These degeneracies leave their footprint directly in the likelihood profiles derived from scattering and oscillation data due to the way the NSI parameters enter into the fundamental observables; therefore, it is important to understand the degeneracy structures in order to know how they can be broken. We outline two such classes of degeneracies between NSI parameters that present themselves in neutrino oscillation and scattering datathose that exhibit degeneracy between NSI parameters of different fermion index f , and those between different flavor indices αβ. Many of the degeneracies we will discuss have already been derived and discussed before in the literature, but we include them here to have a complete motivation of the subsequent analysis.

"Fermion" Degeneracies
The first class of degeneracies concerns the ability for an experiment to distinguish NSI between different SM fermions, f and f , manifested between f,V αβ and f ,V αβ , for example. This type of degeneracy manifests itself differently within three important classes of interactions, namely neutrino oscillations, neutrino-nucleus scattering, and neutrino-electron scattering.

Oscillation Experiments
An experiment measuring neutrino oscillations through the Earth has direct sensitivity to the oscillation Hamiltonian and its NSI contribution to the matter potential; taking 3 flavors α,β = e, µ, τ . The first term, dependent on the neutrino energy E ν , controls flavor oscillations in vacuum. It contains the mixing angles θ 12 , θ 13 , and θ 23 within the PMNS mixing matrix U and the neutrino mass splittings in M = diag[0, ∆m 2 21 , ∆m 2 31 ] which we have taken to be in normal heirarchy (NH). Neutrino NSI are contained in the matter potential V which is a function of the coordinate x; where f,V αβ ≡ f,L αβ + f,R αβ . The matter potential term is responsible for the Wolfenstein effect as well as parametric resonances from neutrino oscillations through the Earth [1,9]. It is dependent on the electron, up-quark and down-quark number densities as well as the NSI parameters. The NSI terms can be expressed as a matrix containing all the flavor αβ vertices; The electron, up, and down NSI enter into the oscillation Hamiltonian in linear combination (Eq. 2.2), and upon diagonalization it is this linear combination that enters into the survival and transition probabilities as physical observables. One may refer to [10] for the treatment of neutrino oscillation amplitudes, and more recently [11] for the exact probability expressions. To simplify things, the electron number density can be factored out and since n u ≈ n d ≈ 3n e in the Earth to good approximation, any experiment that measures the matter potential effects of neutrino oscillation is only sensitive to the sum. We therefore define a phenomenological NSI parameter; This is the NSI observable for experiments measuring NSI in oscillations through the Earth. It defines a plane of solutions in ( e,V αβ , u,V αβ , d,V αβ ) space; unless two out of the three terms are fixed, oscillation experiments have a three-fold degeneracy in their sensitivity to the e, u, and d NSI.

CEνNS Experiments
Neutrinos interact coherently with nuclei if their transferred momentum q satisfies qr n << 1 for a nuclear radius r n . This process is described via the Coherent Elastic Neutrino-Nucleus Scattering (CEνNS) mechanism [12,13]. In this energy range (typically for E ν as high as 100 MeV for most nuclei), the first term in Eq. 2.1 becomes large leaving the matter potential as a subdominant effect and relegating new physics observables to the CEνNS cross-section. The cross-section is given by where we traditionally take the Helm parameterization of the form factor F (q 2 ) with a neutron skin radius r n = 5.5 fm [14,15]. The Q 2 V factor ordinarily contains the SM charges, but with the presence of NSI, potentially allowing flavor changing processes such as ν α + N → ν β + N , it is modified to where s w = sin θ w is the sine of the Weinberg angle and we have taken the initial state flavor α and summed over final state flavors. Z and N are the proton and neutron numbers of the target nucleus, respectively. A linear combination of the up and down vector NSI, which we denote N αβ , can then be factored out such that Q V is a function of a single NSI parameter, which we denote by N αβ ; giving rise to two lines of solutions in the ( u,V αβ , d,V αβ ) plane; The slope of these lines is given by the ratio (2N +Z)/(2Z +N ) which varies depending on the detector material; for instance, 126 54 Xe has a ratio of 1.1 while 40 18 Ar has 1.07. Therefore using CEνNS data from multiple detectors of different materials can have complementary likelihood profiles that can help to break this type of degeneracy [13,16,17].

Elastic Neutrino-Electron Scattering (EνES) Experiments
Lastly, we consider the elastic neutrino-electron scattering (which we will refer to as EνES) crosssection, measured from the Solar neutrino flux, for example, contains a similar charge structure when NSI are included. It may also permit flavor-changing scattering processes ν α + e − → ν β + e − ; where g L = sin 2 θ w − 1 2 and g R = sin 2 θ w . The δ eα δ eβ term encodes the charged-current enhancement for ν α + e − → ν β + e − scattering, and the other Kronecker delta terms take care of removing the SM charges when the process is flavor-changing. There is only one fermion index, e, appearing in this cross-section. In this case we simply identify the phenomenological NSI with the physical NSI; However, the way these NSI appear in the cross-section gives rise to a more complicated degeneracy structure than in CEνNS and oscillation experiments, since the left and right chiral NSI parameters do not simply factor out of the energy-dependent terms in Eq. 2.6. However, we can visualize the degenerate features of this cross-section by equating the NSI cross-section with the SM one and match like-terms in the E r expansion; This system of equations has one or two solutions depending on the presence of the δ eα δ eβ term (and in the two solution case, the larger of these solutions may be disallowed by existing constraints). Note that as m e /E ν → 0, the number of equations reduces to 2 and the number of solutions rises to 4, therefore encouraging the measurement of this cross-section at relatively lower neutrino energies, for example, the 7 Be solar neutrino flux (E ν ≈ 0.86 MeV). For α = β = e only the top line of Eq. 2.7 is a solution, but for all other αβ pairs both solutions exist, implying at most two solutions for e,V αβ . The Eqs. 2.7 are shown in Appendix B for both of the aforementioned cases. We will also examine what happens in the case that more than one NSI flavor index is nonzero at once in the following section. These degeneracy structures have been studied in more detail in the context of the DUNE near detector in [18].
Collecting these three fermion degeneracies together in the simple case of a single flavor index activated, we show the overlapping solutions which are degenerate with the SM in ( u,V ee , d,V ee , e,V ee ) space in Figure 1. These plane solutions correspond to Eqs. 2.3, 2.5, and 2.7. The point at which all solutions simultaneously intersect, i.e., the maximum likelihood point for a likelihood function defined over the combination of oscillation, CEνNS, and EνES data, is at the origin. The important implication here is that the degeracies among triads of ( u,V αβ , d,V αβ , e,V αβ ) NSI parameters can be broken by combining the three aforementioned experimental classes.

"Flavor" Degeneracies
In addition to providing control over degeneracies between e, u, and d NSI, there is another class of degeneracies present to which oscillation and scattering experiments can conspire to resolve, existing between NSI of different flavor vertices, for example, between f,V αβ and f,V γδ .

Oscillation Experiments
and an additional transformation in the off-diagonal NSI parameters, O αβ → −( O αβ ) * , which we do not consider in this work for limiting ourselves to real-valued NSI. These transformations become especially relevant when one also considers the mass-ordering parameters, mixing angles and phases in the vacuum part of the Hamiltonian to vary alongside NSI, giving rise to the generalized mass-ordering degeneracy [6,[19][20][21].

CEνNS Experiments
Allowing more than one NSI of different flavor indices αβ within the Q 2 V factor of the CEνNS crosssection gives rise to generalizations of the plane solutions outlined in the previous section. To see this, once again we equate the NSI-modified Q 2 V with the SM Q 2 V and find all the solutions. For example, taking the incoming neutrino flavor as e, we find the equation

EνES Experiments
Now consider the EνES cross-section but with two NSI flavor indices nonzero, each with left and right chiral components ( e,L ee , e,R ee , e,L eµ , and e,R eµ , for example). Since the current generation of experiments sensitive to EνES are flavor-blind (just like in the CEνNS case), e − + ν e → e − + ν e and e − + ν e → e − + ν µ have indistinguishable final states, in which case we have to sum the cross-section over final states. This gives rise to the following system of equations; One may check that this system of equations has a single real solution e,L ee = e,R ee = e,L eµ = e,R eµ = 0, provided m e /E ν remains of order 1. As m e /E ν tends to 0, the sets of curves defined by the above three equations no longer intersect in NSI space, which translates to a less resolved profile likelihood if one would perform a maximum likelihood estimation on neutrino scattering data sensitive to the EνES cross-section.
What we hope to convey at this stage is that as one continues to introduce more NSI parameters into the scattering cross-sections and oscillation Hamiltonian, the number of solutions may increase and become less resolved in a likelihood analysis, but the relationships between NSI of different fermion indices remains the same. In the spirit of Figure 1, the neutrino oscillation matter potential, the CEνNS cross section, and the EνES cross section form a weak mapping between the phenomenological NSI to physical NSI that exploited by combining the results of the three classes of experiments ( Figure 2). In order to concretely establish sensitivity to the NSI operators that we have considered, it is essential to combine all three of these classes of experiments together in a global analysis.

Combining Oscillation, CEνNS, and EνES Data
Motivated by the degeneracy structures we have just discussed, we will now attempt to illustrate how CEνNS, EνES, and oscillation experiments can be joined together in a global analysis. We will work under the pretense that NSI of all fermion indices f = e, u, d are free to vary; in other words, none of them will be fixed to zero or other values on the basis of external limits. Doing this permits 24 NSI parameters in the count; 12 from e,L αβ and e,R αβ , 6 from u,V αβ and 6 from d,V αβ . As we outlined in Section 1, we take this relatively large set of parameters for two key reasons that we will summarize again. The first reason is that one should not be limited by the pretense that if there is new physics in the neutrino sector, the new vector operators should be restricted to simple combinations of lepton non-universal or flavor-changing neutral currents. The second reason is primarily technological; that we would like to be able to include many NSI parameters (and many experimental data) to understand the full space of degeneracies with computational impunity.
To build our ensemble of experimental data for the global fit, we first begin by considering a minimal setup consisting of a single CEνNS dataset and a single EνES dataset. The solar neutrino spectrum measured by Borexino in Phase II of the experimental program [22] gives us control over the e-NSI through EνES, while the COHERENT collaboration's open data release and observation of the CEνNS interaction at a CsI detector [23] provides sensitivity to the u and d quark NSI. Analyses have already been performed with these datasets [5,24,25], but to explore the full space of NSI degeneracies we will allow an NSI set of maximum size to be free in our analysis, which has not been done before. Secondly, the stopped-pion spallation neutrino source (SNS) at COHERENT only produces ν e , ν µ , and ν µ neutrinos, and with a short baseline to the detectors, there is negligible oscillation of the neutrino flux into τ flavors. This implies a lack of sensitivity to u,V τ τ and d,V τ τ parameters. To supplement our analysis with a dataset sensitive to u,V τ τ and d,V τ τ NSI, we extend the experiments considered so far to include future projections at a future liquid xenon (LXe) dark matter detector (DMD); being optimized for detecting nuclear interactions with the DM halo, a kiloton-scale detector would be sensitive enough to access the CEνNS and EνES cross-sections from naturally occurring neutrino fluxes. This is the "neutrino floor" that DM direct detection experiments are soon set to encounter, and may be enhanced in the presence of NSI [26]. We will project sensitivity to u and d NSI from atmospheric neutrinos, which have the full flavor-range to access u,V τ τ and d,V τ τ NSI parameters through enhancements to the CEνNS cross-section. Such an experiment would also be sensitive to solar neutrinos, which lie at energy ranges such that they mainly interact through neutrino-electron scattering, provided that they can discriminate between electron and nuclear recoils. Therefore we also project the sensitivity to e,L αβ and e,R αβ NSI with solar neutrinos at an LXe DMD to complement the Borexino analysis.
Finally, for the neutrino oscillations aspect of our analysis strategy, we will project simulated data at DUNE from atmospheric neutrinos oscillating through the Earth, which have been studied before in a variety of contexts for measuring the PMNS matrix and NSI [27][28][29]. We stress here that atmospheric neutrinos will supply DUNE with a very rich oscillation dataset to supplement data from beam neutrinos. The atmospheric neutrino flux contains a host of all e, µ, and τ flavor neutrinos after oscillation through the Earth's mantle, and effectively comprises a range of oscillation baselines as neutrinos propagate through the range of zenith angles.

Prior-flow
Now we will outline the statistical treatment for the simulated and measured data in the experiments considered. We take a Bayesian inference approach to constraining the NSI parameters. In the following discussion we use the Bayesian inference package MultiNest [30] to construct likelihood functions and compute the posterior probability distributions of our NSI parameters given the simulated data at each experiment.
To combine the likelihood information from several experiments, traditionally what is done is the construction of a single global likelihood function calculated by the simultaneous simulation of the event spectra for each class of experiment considered. In the context of this analysis, the joint likelihood function would take physical NSI as model inputs, giving L = P (D| , H) for NSI parameters , observed data D = ∪ N i=1 D i for experiments i = 1, . . . , N , and a null hypothesis H for which we take all NSI as zero, i.e. = 0. This style of approach has been taken before in a variety of global analysis settings for neutrino NSI [7,8,31]. This approach can be very computationally expensive, and with 24 parameters in our consideration (12 from e,L αβ and e,R αβ , 6 from u,V αβ , and 6 from u,V αβ ) it may be difficult to accurately discover the posterior distribution in such a large prior volume, let alone converge at all in the evidence computation. This is even without taking into account the potentially many experimental nuisance parameters for detector response, background or signal uncertainties, etc., in the likelihood fit. In fact, in a typical global analysis we need not allow so many NSI parameters to be nonzero at the same time; experimental nuisance parameters can be enough to create a likelihood parameter space of relatively large dimension.
Instead, we take a divide-and-conquer approach illustrated as follows. Suppose we aim to measure NSI parameters x and y. Then suppose we posses or simulate data from two experiments A and B such that experiment A is sensitive to x and y while experiment B is sensitive to y alone. One can then use experiment B to measure y and its posterior distribution given the data at B, π(y | D B ), and subsequently take this posterior distribution as the prior distribution on y for experiment A. In the context of Bayes theorem, where Z is the Bayesian evidence. Since experiment B provides no information on x, we take a uniform prior density u(x) over an appropriate interval such that the joint prior becomes By doing this, we effectively constrain y at experiment A using its prior distribution from experiment B. If A and B have complementarity between any degeneracies that may exist for x and y, they would be combated just as they would by directly calculating π(x, y | D A ∪D B ) via a joint likelihood function for the data from the two experiments A and B. This strategy can be repeated for numerous experiments and with more parameters. By allowing posterior information to "flow" from one experiment into another, we effectively reduce the prior volume that needs to be searched. This is what we propose using data from the COHERENT and Borexino collaborations, in addition to projections for a future LXe DMD and the DUNE experiment, bearing in mind that this scheme could be extended to a variety of others.
The prior ordering structure is shown in Figure 3. At the top level, only uniform priors are used, and by default we fix the uniform interval to (−1, 1) for all vector NSI (and (−0.5, 0.5) for L and R components) for the sake of simplicity. Each subsequent experiment in the "prior-flow" takes its prior from the joint posterior distributions of the experiments above, for the relevant subset of NSI to which those experiments are sensitive, with uniform priors for the NSI that remain. The explicit sets of and their priors are listed in Table 1.

Copulas
With such a strategy, there is an important question as to how one models the prior distributions of a multivariate set of NSI parameters. We remind the reader that according to our strategy, this joint prior of the parameters is constructed based on the posterior distributions of the previous experiments.
If one only uses the one-dimensional marginal distribution as individual prior on each parameter, important correlations between the NSI will be lost. Therefore, we elect to model the joint prior distribution as completely as possible. To do this, we use a copula. In d dimensions, a copula C is a cumulative distribution function (CDF) C : [0, 1] d → [0, 1] with uniform marginal distributions. See [32,33] for a review. Sklar's theorem [34] states that for every d-dimensional joint CDF, in our case F( 1 , . . . , d ) for NSI parameters 1 , . . . , d , there exists a d-copula C such that where F 1 , . . . , F d are the marginal distributions of the NSI parameters. Copula functions, in essence, connect the marginal distributions and the joint distribution through a correlation structure. Given absolutely continuous marginal distributions and the joint distribution, the copula function is unique.
The copula C is usually a function, which can sometimes be written in closed form, whose form is associated with the dependency structures of a known family of statistical distributions. There are many families of copula, and no single copula is guaranteed to be a perfect model of the underlying joint distribution, so in practice one usually chooses the family that best fits the sample data of the joint distribution. For example, the band-shaped degeneracy contours between pairs of u,V αβ and d,V αβ NSI (which we will see in Section 5) may be well-modeled by the Frank family of copulas that capture this kind of correlation well. However, one may also use an empirical copula to fit the joint prior distribution provided one has sample data from MultiNest. We elect to use this option to fit the posterior joint distributions and to subsequently simulate prior distributions in the prior-flow, since this is the most robust and accurate way of modeling the NSI dependency structure discovered by each experiment. To do this we use the R package copula to fit empirical copulas to the joint distributions for each experiment in the prior flow.
Π e Xe = C(F ( e,V αβ )) Figure. 3: The "Prior-Flow" of joint probability information from experiment to experiment. We begin with EνES from solar neutrinos at Borexino and CEνNS at COHERENT, then proceed to EνES and CEνES scattering at a future 1 kton·year LXe dark matter detector, and finally to future atmospheric neutrino oscillation measurements at DUNE. The components of each prior that are inherited from a previous experiment are denoted by Π, which are taken as empirical copulas of the relevant marginals of the NSI parameters. If one of the previous experiments is not sensitive to a particular NSI, the uniform distribution is taken as a prior for that corresponding parameter.
In MultiNest, the prior is formally implemented as a map from the d-dimensional cube [0, 1] d → R d via the inverse CDF or quantile function of the prior distribution. This naturally lends itself to the implementation of copulas, since the available methods of simulation (finding the inverse of a multivariate CDF) are well documented. In order to simulate samples from the empirical copulas, we employ the conditional distribution method [32,35] followed by using the inverse CDFs of the prior marginals to extract sample NSI parameters for each iteration in MultiNest. For an empirical copula C this procedure is illustrated as follows; where the {x i }, i = 1, . . . , d correspond to the d NSI parameters used in each prior and their 1-dimensional marginal CDFs F 1 , . . . , F d . In step 2 we compute the the conditional distributions of the copula, C r|r−1,...,1 , which is equivalent to finding its partial derivatives and taking the pseudoinverse [36]. Since the empirical copula C and its partial derivatives are numerically computed from the MultiNest observations, this inverse is found numerically as well. As an example of how closely the empirical copula can model one of the joint posterior distributions on the NSI, a comparison between the kernel density estimates of the MultiNest samples and the corresponding empirical copula simulation is shown in Figure 4.

Stopped-pion neutrinos at COHERENT
The COHERENT experiment uses the stopped-pion spallation neutrino source (SNS) to produce MeV neutrinos. A CsI detector, situated 19.3 meters away from the source, collects scintillation light to reconstruct the energy E r from CEνNS nuclear recoils in the detector volume.
We use the procedure detailed in Ref. [37] which combines both energy and timing data in the neutrino spectrum at COHERENT from the 4466 kg·days of exposure at the CsI detector. By using probability densities of the time and energy spectra for the ν µ , ν e , andν µ flavor components we can predict the number of observed prompt and delayed neutrino counts. To compute the number of events between recoil energies E a r and E b r , we convolve the neutrino flux for each neutrino species α with the CEνNS cross-section between kinematically allowed neutrino energies; Here we take E min ν = ( E 2 r + 2m N E r + E r )/2 for a nucleus mass m N and neutrino fluxes dΦ να /dE ν for ν µ ,ν µ , and ν e flavor neutrinos.
To take advantage of the multiple detector materials available in by the SNS, in addition the data from the CsI detector we also include a projection of CEνNS data at the liquid argon (LAr) CENNS-10 detector. This exposure is in anticipation of the upcoming CENNS-10 publication, based on the talk New Results from a CEvNS Search with the CENNS-10 Liquid Argon Detector given at Fermilab on January 10, 2020 [38]. In the absence of any LAr data release, we simulate CEνNS data given the null hypothesis (SM) at the LAr detector. Including an argon detector in the analysis will give the CsI CEνNS measurement some complementarity in the ( u,V αβ , d,V αβ ) plane as discussed in Section 2.1. We tentatively take an exposure of 1.5 years (≈ 15.33 ton·days) and an energy range of E r ∈ [20,100] keV in accordance with the expected low-energy threshold of the detector and the upper limit of the coherent scattering regime with the argon nucleus. Energy (or photoelectron) resolutions and time resolutions are taken to be the same as for the CsI detector.
For the likelihood analysis, we again refer to Ref. [37], but instead we take the real-valued u and d NSI as model inputs. We include u,V ee , u,V µµ , u,V eµ , u,V eτ , u,V µτ , d,V ee , d,V µµ , d,V eµ , d,V eτ , and u,V µτ , but not u,V τ τ or d,V τ τ because the negligible presence of τ -flavor neutrinos at the SNS.

Solar Neutrinos at Borexino
The Borexino collaboration has measured the solar neutrino energy spectrum [22] over 92.1 kton·days which provides an important dataset to help constrain NSI in neutrino-electron scattering events. We follow Ref. [39] to model the solar neutrino energy spectrum at Borexino. The solar neutrino event rate is predicted using the EνES cross-section (Eq. 2.6) and convolving it with the oscillated solar neutrino flux; where we assume no direction reconstruction on the incoming neutrino and no flavor-sense, hence, the incoming, oscillated, and transition flavors α β and γ are summed over. The minimum neutrino energy is the same as the one used in Eq. 4.1 but with the replacement m N → m e . We select data about the 7 Be compton edge, corresponding to recoil energies from 550 keV to 1 MeV. As mentioned in Section 2.1, this energy range is sensitive to NSI contributions to the m e /E ν proportional terms in the EνES cross-section, which in turn helps converge on a single NSI solution during the likelihood analysis. This region also contains contributions from radiochemical backgrounds; 85 Kr, 210 Po, 11 C, and 210 Bi.
A log-likelihood function is constructed from the Borexino Phase II data N o i , while the error standard deviation and the expected number of events are denoted by σ i and N s+b i , respectively. The final likelihood combines information from all energy bins i = 1, 2, . . . . We allow the predicted event rate to vary as a function of k = (k Be , k P o , k Kr , k Bi , k C ) which parametrizes uncertainties in the 7 Be flux and background rates. For a background rate R j , we allow it to fluctuate via the parameter k j as follows; We then take the Gaussian prior for these nuisance parameters k with means of 0 and widths given by the rate uncertainties. The predicted event rate is of course a function of NSI as well, taking = ( e,L ee , e,L ee , e,L µµ , e,L τ τ , e,L eµ , e,L eτ , e,R µτ , e,R µµ , e,R τ τ , e,R eµ , e,R eτ , e,R µτ ). The log-likelihood function is now given in Eq. 4.3.  Figure. 5: The solar neutrino spectrum from the Borexino Phase II dataset is shown around the 7 Be compton edge. Data, and information on backgrounds, was obtained from the Borexino data release corresponding to Ref. [22] and Ref [39]. Our standard interactions (SI) prediction is shown in solid red. An example NSI solution that enhances the event spectra is shown in dashed red.

Atmospheric Neutrinos at DUNE
Neutrinos produced in the Earth's atmosphere from cosmic ray processes consist of the ν e ,ν e , ν µ , and ν µ flavor states. These neutrinos are then free to propagate through the Earth and undergo flavor oscillations. The neutrinos can then be detected by the DUNE far detector, capable of reconstructing the neutrino energy and direction (or zenith angle between the incoming neutrino trajectory and the horizon plane at the detector).
Since we limit the scope of this analysis to neutral-current vector NSI, we describe the chargedcurrent interactions in LAr with the SM prediction, in particular ν α + n → − α + p + ν α + p + → + α + n via charged-current quasi-elastic (CCQE) scattering. We ignore resonance production processes, which will reduce statistics but mitigates the theoretical uncertainties in the resonance production crosssections as well as the hadronic energy corrections that smear the energy reconstruction [40]. Additionally, we will restrict ourselves to ν µ scattering, whose final state typically gives rise to two well-identified charged tracks in the detector. For the analytic form of the cross-section σ(E ν ) we implement the one developed by in Ref. [41] which includes a parameterization of the transverse enhancement from meson exchange currents inside the nucleus. This parameterization offers a decent fit to cross-section data at in the relevant energy range for atmospheric neutrinos of 100 MeV to 1 GeV. We provide more details of the implementation of this cross-section in Appendix A.
Atmospheric fluxes Φ α (cos θ, E ν ) are taken from the FLUKA results of Ref. [42] for the Super-Kamiokande site. To obtain the predicted event count of a neutrino flavor α between energies E a ν and E b ν and zeniths cos θ 1 and cos θ 2 , assuming perfect reconstruction, we convolute σ(E ν ) with the oscillated atmospheric flux. The number of neutrinos observed for a flavor α is then given in Eq. 4.4:

.4)
P βα = P βα (cos θ, E ν ) are the survival and transition probabilities determined from the neutrino Hamiltonian with NSI. To calculate P βα we employ a numerical diagonalization method on the oscillation hamiltonian in conjunction with a simplified version of the PREM model for the electron number density in the Earth. The rich spectrum is shown in Figure 6 for all zenith angles and energies between 100 MeV and 1 GeV. The most NSI-sensitive region lies below the horizon, so we select 20 zenith bins for cos θ ∈ [−0.975, −0.025]. This corresponds to an angular resolution of about 18 • . We use 20 energy bins between 100 MeV and 1 GeV. We take a 10 year exposure with the full 40 kton far detector volume. We then employ the same method as in Eq. 4.3 but now defining the log-likelihood function over both energy and zenith bins. Since neutrino oscillations are sensitive to e,V αβ , d,V αβ , and d,V αβ in the Earth's matter potential, we allow all 18 real NSI degrees of freedom to vary in the likelihood function.

Solar and Atmospheric Neutrinos at a Future LXe Dark Matter Detector
The atmospheric and solar neutrino event rate at a LXe dark matter detector via CEνNS (Eq. 2.4) or EνES (Eq. 2.6), respectively, can be predicted via a similar convolution to Eq. 4.4. For the LXe detector we assume no direction reconstruction on the incoming neutrino and no flavor-sense; just energy reconstruction via nuclear recoils. Therefore we only use the zenith-integrated flux and sum over incoming neutrino flavors. For the statistical analysis of the predicted data, we again use a loglikelihood as in Eq. 4.3. Once again, for solar neutrinos we use a sum over energy bins and we include 12 NSI degrees of freedom (6 e,L αβ and 6 e,R αβ ) just as in the Borexino analysis. Since atmospheric neutrinos may oscillate into τ flavors in the Earth and interact in the detector via CEνNS, we are now sensitive to u,V τ τ and d,V τ τ , thereby expanding our NSI degrees of freedom from the set used in the COHERENT analysis from 10 to 12 NSI.
We set the design goal exposure for this future detector to be 1 kton·year. While this is larger than existing proposals in the literature [43], we take the approach of understanding the physics reach of such an experiment for an optimistic exposure. For the energy threshold, we assume recoils can be realistically reconstructed as low as 5 keV, looking for CEνNS events up to 50 keV. To see CEνNS events from atmospheric neutrinos, again we take the FLUKA result for the atmospheric flux, but an important point needs to be raised regarding this flux; atmospheric neutrinos need to have low enough energies ( 60 MeV) to scatter coherently off Xe nuclei, but the 3D FLUKA result only goes as low as 106 MeV. Therefore we extrapolate the atmospheric neutrino FLUKA fluxes down to 10 MeV using a 3rd-order polynomial in log space, shown in Figure 7. While calculations of the zenith-by-zenith flux do not exist yet down to 10 MeV, we can check that the zenith-integrated flux agrees well with the one reported in in Ref. [44] for the solar-averaged flux at Super-Kamiokande. We use a 3rd-order spline fit in log space in order to perform the extrapolation, i.e., log Φν = α + β log Eν + γ(log Eν ) 2 + δ(log Eν ) 3 for fit constants α, β, γ, δ.
There is one final remark; this class of detector would also be sensitive to CEνNS from solar neutrinos coming from the 8 B processes, inducing nuclear recoils up to energies of a few keV. For the analysis in this paper, we do not include this contribution to the event rate, in order to focus on atmospheric neutrino-induced events greater than 5 keV energy recoil. Extracting the recoils from 8 B neutrinos would require a more dedicated analysis, as the complete detector efficiencies are difficult to estimate at this stage.
A quantitative summary of the specifications for each experiment simulated in this section can be found in Table 2.

Results
Using the methods we have just described, we derive the fits to the NSI parameters at each stage of the "prior-flow" outlined in Figure 3, with a final joint posterior distribution derived from the last stage at DUNE. We show the 1-dimensional marginalized posterior distributions for the 18 real-valued vector NSI parameters at various stages of the prior-flow in Figure 8. Good convergence on the electron NSI e,V αβ is observed, but we note that the posterior means for e,V ee and e,V µµ are slightly negative to accommodate the best-fit on the Borexino data about the 7 Be edge. We state the 95% credible intervals corresponding to the 18 electron, u and d quark NSI in Table 3. In the table we compare the credible intervals for COHERENT and Borexino (middle column) with those for the projections at DUNE and the LXe DMD. We observe that DUNE and the LXe DMD make an improved reduction in the width of the credible intervals on electron NSI by a factor of 2 to 3.
Convergence for u,V αβ and d,V αβ is also improved by DUNE and the LXe DMD with respect to the posteriors from COHERENT. COHERENT and the LXe DMD offer good constraints on the u and d quark NSI that enter in as priors for DUNE, but the constraints on electron NSI from solar neutrinos also help constrain the u and d quark NSI indirectly via the linear correlation E ee = e,V αβ +3 u,V αβ +3 d,V αβ , which enters into the matter potential to which DUNE is sensitive. Phenomenologically speaking, a strong constraint on e,V αβ incurs an equal and opposite constraint on 3 u,V αβ + 3 d,V αβ . Note, however, that this relationship can also have the effect of inducing biases; if the data at DUNE is consistent with E ee = 0, then via the aforementioned linear combination, any bias in e,V ee will induce a bias in u,V αβ and d,V αβ via their correlations through E ee . We see this effect notably in ee and µµ NSI caused by the negative bias in e,V ee from the Borexino part of the analysis. It should also be pointed out here that the biases seen in Figure 8 are reflected in the credible intervals in Table 3; for some NSI, the fit has pushed the credible interval to exclude the zero value point, but again this arises as an artifact of the null hypotheses we have assumed for DUNE and the LXe DMD and the intrinsic correlation between the fits on e,V ee , u,V αβ , and d,V αβ NSI. Additionally, even with multiple detector materials available in our analysis to break the u,V αβd,V αβ degeneracy, some degeneracy still remains from the correlation between u and d quark NSI in the 2-dimensional marginal posterior distributions; see Figure 9 where we show all the 2-dimensional projections of the prior-flow posteriors in 18 NSI dimensions. The credible regions in the ( u,V αβ , d,V αβ ) planes are certainly improving with each stage, but the correlation never fully goes away. Notice also that in the flavor-diagonal u and d quark NSI in Figure 9, single degeneracy bands are seen for the COHERENT analysis (yellow), while we see a double-band degeneracy for LXe in the ( u,V τ τ , d,V τ τ ) plane. This is because in the COHERENT CsI data, a deficit in CEνNS events was observed, leading to a modification of Eq. 2.5 that brings the second line of solutions closer to the line intersecting (0, 0). On the other hand, since we simulated data at the LXe DMD based on the SM expectation for CEνNS, there is no deficit in the event rate and therefore the two solution lines in accordance with SM degeneracy are in full view for ( u,V τ τ , d,V τ τ ). This double-solution degeneracy is nicely broken by DUNE which can be seen in Figure 9.
By defining q,V αβ ≡ u,V αβ + d,V αβ we can transform away the strong correlation between u and d quark NSI and visualize the remaining degeneracy between electron and quark NSI that DUNE would exhibit. In Figure 10 we plot a grid of the 1-and 2-dimensional marginal projections of the NSI parameters reduced to just 12 NSI (6 e,V αβ and 6 q,V αβ ). Again the double-solution degeneracy on q,V τ τ is broken by DUNE. After this transformation we observe good convergence on q,V αβ and e,V αβ , with improved reduction in the credible interval widths by roughly a factor of 2 with the addition of DUNE and the LXe DMD as shown in Table 3. Some correlation remains between the pairwise combinations of q,V αβ and e,V αβ , most prominantly between q,V µτ and e,V µτ . Overall, the reduced set of 12 NSI comprising e,V αβ and q,V αβ has the best convergence with the most number of degeneracies broken, while still representing a set of NSI parameters that are not too phenomenological to be non-influential to model-building.
Finally, we also show in Figure 11 the posteriors for e,L αβ and e,R αβ NSI before they are passed in as priors for DUNE in their vector combinations. The 68% credible contours and 1-dimensional marginal posterior probability distributions are compared between Borexino and a future LXe DMD. Excellent convergence is achieved on e,L ee and e,R ee due to the CC enhancement to the EνES cross-section which constructively interferes with e,L ee and e,R ee NSI to produce larger effects on the 7 Be flux.    αβ , e,R αβ ) for a total of 12 NSI degrees of freedom. The distributions and contours for Borexino (yellow) and a future LXe DMD with priors from Borexino (magenta) are overlayed.

Conclusion
We have shown that it is possible to measure neutrino NSI, significantly breaking their degeneracies, even when many NSI parameters are nonzero. The inclusion of three different classes of observables -the CEνNS and EνES processes and neutrino oscillations -are essential to constructing a global analysis whose experimental data are complementary to one another in the NSI model parameter space. We have chosen COHERENT and Borexino datasets as excellent representative neutrino scattering datasets, but these can readily be augmented with a variety of others. The far detector at DUNE with its large volume should provide excellent constraints on NSI through its ability to access rich oscillation information through the detection of atmospheric neutrinos after they interact with the matter potential of the Earth. The addition to this ensemble of neutrino scattering data at future dark matter experiments we showed to be a natural complement to the CEνNS data at COHERENT by their potential sensitivity to τ flavor neutrinos from solar and atmospheric sources. We stress that the experiments considered here are best used together as a unified source of data to investigate neutrino NSI.
The relatively many NSI considered in the analysis and multiple experiments being simulated became pragmatically realizable with our divide-and-conquer approach using the copula. We demonstrated that a strategy of connecting posterior probability distributions as Bayesian priors from experiment to experiment allows one to scale a global analysis with a potentially large number of model and nuisance parameters, with copulas facilitating the transfer of prior information. This novel "priorflow" framework we outlined can be extended in a straightforward way to include other existing data which would be sensitive to NSI. The Bayesian estimation of posterior probability distributions on the relatively large number of NSI parameters considered here was demonstrated to be tractable.
Our analysis could be extended, notably, to include neutrino-nucleus scattering data from CHARM, whose measurement of the cross-section ratio of NC to CC processes provides a well-known complementary constraint to CEνNS measurements in the ( u,V αβ , d,V αβ ) plane. It was omitted from this work for not providing a strong enough constraint relative to the parameter ranges we restricted ourselves to ( f,V αβ ∈ [−1, 1]), but for a broader parameter space it would be interesting to integrate CHARM data into the analysis strategy [20]. Additionally, there are numerous oscillation data sets readily available which could contribute to the statistical power of the analysis, integrated into a global NSI study in a similar manner to DUNE.
To further generalize the projected constraints on NSI, we plan to extend our investigation to include complex-valued NSI parameters as well as effective NSI operators in scenarios where the underlying mediator masses are light and comparable to the scale of the neutrino momentum transfer. We emphasize the importance to obtain constraints on the NSI in these more general scenarios in order to support model-independent results and drive more theoretical work in this area. Since the space of neutrino experiments is expanding quickly and allowing for highly comprehensive analyses in the future, the need for new tools to combat model parameter degeneracies in highly generalized settings will be highly sought after. It is precisely these degeneracies that should make the reader appreciate that neutrino scattering and oscillation experiments should be thought of together, as a unified source of experimental information on new physics. We hope to have cut a pathway with the unique strategy presented here to give global analyses of neutrino interactions the ability scale up as we enter the precision frontier of neutrino physics.

A CCQE Cross-section
We will now review the charged current quasi-elastic (CCQE) scattering process in liquid 40 18 Ar. This concerns the reactions ν α + n → − α + p + andν α + p + → + α + n taking place with the protons and neutrons in the nucleus. For the cross-section and form factors we refer to Refs. [45,46]. Ref. [47] also provides a very comprehensive review but note that equation 57 has incorrectly flipped the sign assignment for ν andν scattering cases.
The CCQE differential cross-section as a function of the momentum transfer Q is where the + sign is taken forν and − for ν scattering. E ν is the energy of the initial state neutrino, M is the target nucleon mass, and s − u = 4M E ν − Q 2 − m 2 where m is the mass of the final state lepton. Let τ ≡ Q 2 /4M 2 and let ξ ≡ µ p − µ n = 4.706, with m being the outgoing lepton mass and M the target nucleon mass (neutron or proton). We take the global fit to the axial mass M A = 1014 MeV. The A, B, and C terms are given as follows Each of the form factors above can be constructed from dipole terms; and we will take G E = G D and G M = ξG D . The transverse enhancement from meson exchange currents [45] is parameterized by where the best fit parameters are a = 6 × 10 −6 MeV −2 , b = 3.5 × 10 5 MeV 2 . The form factors can now be given.
Finally, we scale the total cross-section by the number of target nucleons in the 40 18 Ar nucleus; by Z = 18 forν and by N = 22 for ν scattering. The cross-section per nucleon is plotted for each neutrino type and compared with NOMAD [48] and MiniBooNE [49] data in Figure 12.  Figure. 12: The charged current quasi-elastic (CCQE) cross-section per nucleon is plotted by integrating Eq. A.1 for each neutrino species. Only the νµ scattering cross-section is used in this work to predict νµ scattering rates at DUNE, for energies between 100 and 1000 MeV. The NOMAD [48] and MiniBooNE [49] measurements of the νµ cross-section are overlayed.

B More on EνES Degeneracies
In Figure 13 a visualization of the degeneracy structure between e,L αβ and e,R αβ is shown. The curves shown are defined by equating the constant term (blue), the E r terms (green), and the E 2 r terms (red) in the EνES cross-section with their SM forms. In the case that the initial state neutrino is of electron flavor, the charged current enhancement leads to the solid blue circle which only intersects the other two curves once at the black point at the origin. For other initial state flavors, we obtain the dashed blue circle which intersects the other curves at both black points -a two-solution degeneracy with the SM.