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

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 JHEP09(2020)106 dimension-6 neutral-current (NC) vector NSI among the first-generation of SM fermions with real-valued couplings. 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 modelindependent 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. Although we will not consider complex-valued NSI, the phases φ f,P αβ in the decomposition f,P αβ = | f,P αβ |e iφ f,P αβ exhibit more pronounced degeneracies with the PMNS parameters such as the CP-violating phase δ CP and other vacuum parameters [8,9], which we consider fixed in this work. We will restrict our focus to degeneracies between real NSI and reserve a maximally-general treatment of complex-valued parameters and their degeneracies for a future analysis.

To address issue (ii), we have developed a new statistical technique by virtue of divideand-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.

JHEP09(2020)106

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 data -those 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.

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 (∆m 2 31 > ∆m 2 21 ). 1 Neutrino NSI are contained in the matter potential V which is a function of the coordinate x;

where f,V αβ ≡ f,L αβ + f,R αβ . This potential term supports Wolfenstein oscillations in matter [1], as well as richer phenomena, such as the Mikheyev-Smirnov-Wolfenstein (MSW) effect [10,11] and parametric resonances of neutrino oscillations through discretely varying matter densities [12], e.g., the Earth. 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 [13] for the treatment of neutrino oscillation probabilities, and more recently [14] for the exact amplitudes. 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. 1 In this work we fix θ12 = 0.576, θ13 = 0.148, θ23 = 0.722, ∆m

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 [15,16]. 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 [17,18]. The Q 2 V factor ordinarily contains the SM charges, but with the presence of NSI, potentially allowing flavor changing processes such as

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 us our second phenomenological NSI parameter. This time the linear combination of up and down pieces is responsible for degeneracy between u,V αβ and d,V αβ NSI in CEνNS experiments. Upon squaring Q V and setting N αβ equal to a constant reveals a set of solutions in the ( u,V αβ , d,V αβ ) plane. For example, consider a single flavor-diagonal

to find the solutions degenerate with the SM gives

giving rise to two lines of solutions in the (

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 [16,19,20].

Elastic Neutrino-Electron Scattering (EνES) experiments

Lastly, we consider the elastic neutrino-electron scattering (which we will refer to as EνES) cross-section, 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. Unlike the CEνNS differential cross-section, there is only one fermion index appearing in eq. (2.10) and since e,R αβ and e,L αβ appear next to terms of different energy dependence, it is in principle possible to disentangle them in a multiparameter fit to data. For this reason, 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.10). We can visualize the degenerate features of this cross-section by equating the NSI crosssection with the SM one and match like-terms in the E r expansion;

(2.12) 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). 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 [21].

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.4), (2.9), and (2.13). 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.

CEνNS experiments

Allowing more than one NSI of different flavor indices αβ within the Q 2

V factor of the CEνNS cross-section 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

Expanding eq. (2.17) into the u and d components yields the equation of a hyperellipse in 6 dimensions, but depending on the number of NSI included, this hyperellipse breaks down into simpler solution sets. For example, in the case that we have u,V ee , d,V ee , and u,V eµ nonzero, eq. (2.17) becomes the equation for an infinite elliptic cylinder in three dimensions.

If we add in d,V eµ , cancellations through the linear combination u,V eµ

eµ become available, once again yielding sets of plane solutions;

JHEP09( 2020)106

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 JHEP09(2020)106 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 αβ . It is reminded that we have limited these NSI to be real-valued, and the inclusion of complex phases will be reserved for a later study which also involves the vacuum parameters such as the PMNS phase δ CP . 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 [27] 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 [28] provides sensitivity to the u and d quark NSI. Analyses have been performed with these data sets in both minimal NSI scenarios and for a broad range of operators [5,[29][30][31][32][33][34][35], but to explore the full space of NSI degeneracies we will allow an NSI set of maximum size to be free to vary, 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 explore NSI through the CEνNS and EνES cross-sections from naturally occurring neutrino fluxes [36][37][38][39]. This is the "neutrino floor" that DM direct detection experiments are soon set to encounter, and may be enhanced in the presence of NSI [40]. 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 [41,42]. We stress here that atmospheric neutrinos can supply DUNE with very rich JHEP09(2020)106 oscillation data below E ν = 1 GeV [43] 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. In addition, although it is not considered in this work, the hierarchy degeneracies and the interplay of the CP phase δ CP in the PMNS matrix with NSI can be explored in DUNE [23].

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 [44] 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,45]. 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 1. NSI parameters used in this analysis and the prior scheme used for each experiment.

Here U N are N -dimensional uniform priors on the NSI vector of length N , chosen to range from -0.5 to 0.5. Π e B , Π u,d C , and Π u,d Xe are the priors taken from the posterior distributions at Borexino, COHERENT, and a future LXe DMD, respectively. 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 [46,47] for a review. Sklar's theorem [48] 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

Π 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.

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 captures 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.

JHEP09(2020)106

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 [46,49] 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 conditional distributions of the copula, C r|r-1,...,1 , which is equivalent to finding its partial derivatives and taking the pseudo-inverse [50]. 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 MultiNest samples and the corresponding empirical copula simulation is shown in figure 4.

Stopped-pion neutrinos at COHERENT

There have been numerous analyses of COHERENT neutrino data for the determination of NSI, neutrino vacuum parameters and nuclear structure [5,17,29,32,34,35,38,51,52]. Particularly, we use the procedure detailed in ref. [53] 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 density functions (PDFs) 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 and detector efficiency η(E r );

Here we take Posterior kernel density contours and their marginals plotted for the e,V ee -e,V µµ NSI. The joint distribution observations from MultiNest (top-left) are compared with 10 4 simulated samples drawn from fits for a Gaussian copula (top-right), a Gumbel copula (bottom-left), and an empirical copula (bottom-right). In each fit, the 1-D marginals are reproduced well, but only the empirical copula can accurately reproduce the non-trivial dependency structure of the MultiNest observations in this example.

To take advantage of the multiple detector materials available in by the SNS, in addition the data from the CsI detector we also include data from the liquid argon (LAr) CENNS-10 detector [54]. 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. Spectra for both the CsI and LAr analysis are shown in figure 5.

In the CENNS-10 data release [55], data is binned not only in recoil energies E r and trigger times t, but also in a third dimension, F 90 , corresponding to the light yield fraction in the first 90 ns of the PMT response. Therefore, we predict the 3D binned events in LAr for the using both timing and F 90 PDFs convoluted together with the energy response as in eq. (4.1).

For the CsI likelihood analysis, we again refer to ref. [53] and use a binned log-likelihood that is marginalized over nuissance parameters for the backgrounds (steady-state (SS) and beam-related neutrons (BRN)) and systematic uncertainties.

For the LAr likelihood analysis, alternate PDFs are provided in the data release to encapsulate systematic uncertainties in the BRN and CEνNS rates. This allows us to add parameters to vary the expected CEνNS and BRN rates on top of the NSI variation as follows. First, there are several normalization uncertainties; 13% uncertainty on the CEνNS normalization, 30% and 100% uncertainty on the prompt and delayed BRN components, respectively, and a modified statistical uncertainty of N SS /5 for the steady-state (SS) JHEP09(2020)106 background (in accordance with measurements on SS done in a ×5 oversampled time window). There are several systematics uncertainties; in the BRN timing mean and width, in the BRN energy distribution, in the CEνNS timing mean, and the CEνNS F 90 -Energy distribution. For each of these, there is an alternate PDF template that represents the ±1σ deviation from the BRN or SM CEνNS expectation. We represent the counts in these 3D PDFs over E r , F 90 , and t as vector n over the flattened bins i = 1, . . . , M . Then, for each systematic uncertainty, we parameterize fluctuations in the expected PDFs as transformations of bin contents according to a normal (skew-normal) distribution given by the symmetric (asymmetric) systematic deviations in each bin;foot_0

where F -1 is the inverse-CDF of the normal wit mean 0 and standad deviation of σ i (or skew-normal distribution with the location parameter 0, a scale of σ i and the skewness parameter α i ), whose argument u ∼ U (0, 1) is a uniform variate allowed to float in the MultiNest likelihood scan. After the likelihood scan, the systematic and statistical nuisance parameters will be marginalized over and the NSI posterior distributions will be passed in as priors into the LXe fit as described in figure 3.

For the NSI parameters in both likelihoods, 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 [27] over 92.1 kton•days which provides an important dataset to help constrain NSI in neutrino-electron JHEP09(2020)106

scattering events. We follow refs. [33,56] to model the solar neutrino energy spectrum at Borexino. The solar neutrino event rate is predicted using the EνES cross-section (eq. (2.10)) and convolving it with the oscillated solar neutrino flux;

where we assume no direction reconstruction on the incoming neutrino and no flavorsense, 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.5).

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). In particular, in ref. [43] the atmospheric neutrino flux below E ν = 1 GeV has been shown to exhibit rich oscillation information, not only aids in the measurement of the leptonic mixing angles and the CP violating phase δ CP , but also in gaining sensitivity to the NSI matter effects. 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. [27] and ref. [56]. Our standard interactions (SI) prediction is shown in solid red. An example NSI solution that enhances the event spectra is shown in dashed red.

Since we limit the scope of this analysis to neutral-current vector NSI, we describe the charged-current interactions in LAr with the SM prediction, namely the processes

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 cross-sections as well as the hadronic energy corrections that smear the energy reconstruction [57]. 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. [58] which includes a parameterization of the transverse enhancement from meson exchange currents inside the nucleus. This parameterization offers a decent fit to crosssection 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. [59] 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.6): The rich spectrum is shown in figure 7 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.5) 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.5)) or EνES (eq. (2.10)), respectively, can be predicted via a similar convolution to eq. (4.6). 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 log-likelihood as in eq. (4.5). 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 [60], 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 JHEP09(2020)106 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 8. 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 ref. [61] for the solar-averaged flux at Super-Kamiokande.

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.

JHEP09(2020)106 5 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 9. In addition, the 95% credible intervals corresponding to the 18 electron, u and d quark NSI are listed in table 3. 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. 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 9 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 10 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.

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 11 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 αβ , JHEP09(2020)106 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 12 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.

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 data sets as excellent representative neutrino scattering data sets, 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 JHEP09(2020)106 q,V ee a potentially large number of model and nuisance parameters, with copulas facilitating the transfer of prior information. This novel "prior-flow" 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 αβ )

e,L ee αβ , 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. 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 [24]. 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 JHEP09(2020)106 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. Figure 13. 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 [65] and MiniBooNE [66] measurements of the ν µ cross-section are overlayed.