Improved global fit to Non-Standard neutrino Interactions using COHERENT energy and timing data

We perform a global fit to neutrino oscillation and coherent neutrino-nucleus scattering data, using both timing and energy information from the COHERENT experiment. The results are used to set model-independent bounds on four-fermion effective operators inducing non-standard neutral-current neutrino interactions. We quantify the allowed ranges for their Wilson coefficients, as well as the status of the LMA-D solution, for a wide class of new physics models with arbitrary ratios between the strength of the operators involving up and down quarks. Our results are presented for the COHERENT experiment alone, as well as in combination with the global data from oscillation experiments. We also quantify the dependence of our results for COHERENT with respect to the choice of quenching factor, nuclear form factor, and the treatment of the backgrounds.


Introduction
Experiments measuring the flavor composition of neutrinos produced in the Sun, in the Earth's atmosphere, in nuclear reactors and in particle accelerators have established that lepton flavor is not conserved in neutrino propagation, but it oscillates with a wavelength which depends on distance and energy. This demonstrates beyond doubt that neutrinos are massive and that the mass states are non-trivial admixtures of flavor states [1,2], see ref. [3] for an overview.
Under the assumption that the Standard Model (SM) is the low-energy effective model of a complete high-energy theory, a completely model-independent parametrization of the possible effects of New Physics (NP) at low energies is through the addition to the SM Lagrangian of higher-dimensional operators which respect the SM gauge group. The only dimension-five operator that can be built using just SM fields is the Weinberg operator [4], which coincidentally gives rise to neutrino masses. Given that the observation of neutrino masses is one of the strongest indications of physics beyond the SM, one may therefore expect additional effects from higher dimensional operators. At dimension six, the allowed set is larger and includes four-fermion operators affecting neutrino production, propagation and detection processes, usually referred to as Non-Standard neutrino Interactions (NSI). For example, the effective Lagrangian 2G F ε f f ,P αβ (ν α γ µ P L β )(f γ µ P f ) + h.c. (1.1)

JHEP02(2020)023
would induce non-standard charged-current (CC) production and detection mechanisms for neutrinos of flavor α, while 2G F ε f,P αβ (ν α γ µ P L ν β )(f γ µ P f ) + h.c. (1.2) would lead to new neutral-current (NC) interactions with the rest of the SM fermions. In both eqs. (1.1) and (1.2), f and f refer to SM fermions, denotes a SM charged lepton and P can be either a left-handed or a right-handed projection operator (P L or P R , respectively). Note that the new interactions may induce lepton flavor-changing processes (if α = β), or may lead to a modified interaction rate with respect to the SM result (if α = β). While CC NSI are severely constrained from the precise measurement of meson and muon decays (see for example refs. [5][6][7]), constraining NC NSI is a much more daunting task due to the technical challenges in the computation of neutrino-nucleus interactions, and to the experimental difficulties related to the measurement of neutrino NC interactions. In this case, one may expect to see an observable effect in neutrino oscillations, without entering in conflict with other experimental constraints. Of course, if the set of d = 6 operators is obtained from a NP model at high energies, electroweak gauge invariance generically implies that the NC NSI operators can only be generated together with similar operators involving charged leptons, for which the experimental constraints are much tighter [8,9]. However, recently it has been argued that viable NP models with light mediators (i.e., below the electroweak scale) may lead to large NC NSI effects which would affect neutrino oscillation experiments without spoiling the precise determination of charged lepton observables [10][11][12][13][14][15]. For a recent review on viable NSI models from light mediators see, e.g., ref. [16]. 1 In this case, the best model-independent constraints available in the literature for the vector operators inducing NC NSI come from global fits to oscillation data, which are very sensitive to modifications in the effective matter potential [18,19] felt by neutrinos as they propagate in a medium. Since such modifications arise from a coherent effect, oscillation bounds apply even to NSI induced by an ultra light mediators, as long as their mass is M med 1/R Earth ∼ O(10 −12 ) eV [20]. In particular, oscillation experiments are sensitive to the combinations of Wilson coefficients that is, to vector NSI. It should be noted that, while oscillation data are sensitive to all flavor-changing NSI, only the differences between flavor-diagonal NSI parameters induce observable changes in the matter potential. Consequently, oscillation experiments can only bound five combinations of vector NSI: two diagonal ε f αα − ε f ββ , and three non-diagonal ε f αβ with α = β. Furthermore, in the presence of NSI [18,21,22] a degeneracy exists in oscillation data, leading to a qualitative change of the lepton mixing pattern. This was first JHEP02(2020)023 observed in the context of solar neutrinos, where for suitable NSI the data can be explained by a mixing angle θ 12 in the second octant, the so-called LMA-Dark (LMA-D) [23] solution. This is in sharp contrast to the established standard MSW solution [18,19], which requires a mixing angle θ 12 in the first octant.
The origin of the LMA-D solution is a degeneracy in the oscillation probabilities due to a symmetry of the Hamiltonian describing neutrino evolution in the presence of NSI [24][25][26][27]. For neutrino oscillations in vacuum, it can be easily shown that a simultaneous change in the neutrino mass ordering (that is, the sign of ∆m 2 31 ≡ m 2 3 −m 2 1 , with m i being the masses of the three neutrino mass eigenstates), the octant of the solar mixing angle θ 12 and a shift in the leptonic CP-violating phase δ CP leaves the oscillatory pattern of neutrinos completely unaffected [26,27]. Although this degeneracy is broken in presence of a standard matter potential, this is no longer the case if NSI are allowed [24,25,27]. Hence, the LMA-D degeneracy makes it impossible to determine the neutrino mass ordering by oscillation experiments alone [27], and therefore jeopardizes one of the main goals of the upcoming neutrino oscillation program.
For oscillation experiments performed in matter with a uniform neutron/proton ratio the LMA-D degeneracy stands as exact. Possible ways to lift it are through the combination of oscillation data obtained in environments with different chemical compositions, or for experiments observing neutrino oscillations in a matter potential whose neutron/proton ratio changes sizeably along the neutrino path (as in the case neutrino propagation inside the Sun). This also implies that the results depend on the relative strength of the couplings to up and down quarks. However, global fits to present data show that the LMA-D solution is still pretty much allowed by oscillation data alone [25,28].
As discussed in refs. [23,27,[29][30][31], non-oscillation data (such as those collected by neutrino scattering experiments) is needed to break this degeneracy. However, for light enough mediators (M med O(10) GeV) bounds from deep inelastic neutrino scattering experiments, such as CHARM [32] and NuTeV [33], can be successfully avoided [12,27,30]. Thus, if NSI is generated by mediators as light as about 10 MeV the degeneracy can only be broken through the combination with results on coherent neutrino-nucleus scattering [34] (CEνNS). In fact, CEνNS has been recently observed for the first time by the COHERENT experiment [35] using neutrinos produced at the Spallation Neutron Source (SNS) sited at the Oak Ridge National Laboratory. The analysis of COHERENT data allows to constrain two of the three flavor-diagonal NSI operators, since the neutrino flux contains both muon and electron neutrinos. In combination with the oscillation analysis, this allows for the independent determination of all the three flavor-diagonal coefficients and has a strong impact on the LMA-D degeneracy.
To this end, in ref. [29] we combined the first results from the COHERENT experiment [35] with previous bounds on NSI from a global fit to oscillation data [25]. The outcome of such analysis showed that it was already possible to reject the LMA-D degeneracy beyond 3σ in models where the NSI operators involved only a single quark flavor (either up or down). In addition, the study presented in ref. [29] also yielded the first independent determination of the three flavor diagonal NSI couplings to up or down quarks. Subsequently, in ref. [28] we updated our analysis of oscillation experiments to account for JHEP02(2020)023 the most recent data which had become available, and we also generalized our framework to account for a NSI operators with both up and down quarks at the same time, under the restriction that the neutrino flavor structure of the NSI interactions is independent of the quark type. More explicitly, in ref. [28] we assumed that the Wilson coefficients for the vector couplings to quarks could be parametrized in terms of a quark-indepenent matrix ε η αβ and an angle η as The analysis in ref. [28] concluded that the LMA-D degeneracy is a common feature of the results of the oscillation analysis for a wide range of values of η, i.e., for a broad spectrum of up-to-down NSI coupling strengths. It also showed that the COHERENT data from ref.
[35] was able to lift the degeneracy for a large subset of NSI models. More recently the COHERENT collaboration released the detailed timing and energy information [36] of their data. This allows for further tests of NP and it has been used to this effect in a series of recent works [37][38][39][40][41][42][43][44][45][46][47]. With this motivation, in this work we quantify the effect of including the timing and energy information of COHERENT data on the bounds on NC NSI and, in particular, on the global analysis in combination with oscillation data. In doing so we consider the following novel aspects in the analysis of COHERENT energy and timing data with respect to the analyses performed by other groups [37][38][39][40][41][42][43][44][45][46][47]: • A broader range of NSI models, with couplings parametrized in eq. (1.4), has been considered.
• Following recent discussions in the community regarding the validity of the quenching factor (QF) assumed by the experimental collaboration, in our simulation of the COHERENT signal we use a variety of QF implementations. In particular we introduce a new parametrization from our own fit to the calibration measurements performed by the TUNL group [35,36], never considered before in the literature. The impact of the QF parametrization on the fit is therefore quantified and clarified.
• We study the effect of the nuclear form factor employed and in doing so we use, for the first time, the results of a state-of-the art theoretical calculation of the nuclear form factor provided to us by the authors of refs. [48,49].
• Quantitatively, the most relevant novelty in our analysis is our reevaluation of the steady-state background, as opposed to the ad-hoc parametrization used in the experimental data release (and employed in one way or another in the subsequent phenomenological analyses by different groups). As we will show, using that parametrization results into the disfavouring of the Standard Model at the 2σ level. We show how this is not the case when we use their own beam-off data to determine their steady-state background.
In addition, as mentioned above, we present for the first time the results of combining the information from the COHERENT time and energy spectrum with that from the analysis JHEP02(2020)023 of oscillation data. We quantify the improvement of the bounds on the NSI coefficients by comparing with the results of ref. [28]. Our results show that the LMA-D solution is more significantly disfavoured for a wider range of NSI models. Furthermore, our scrutiny of the possible variations of the COHERENT time and energy information (described above), puts the final bounds on more solid ground. This paper is organized as follows. We start by briefly summarizing in section 2 the framework for the evaluation of the CEνNS event predictions at the SNS in the presence of NSI. Section 3 describes in detail our implementation of the predicted timing and energy dependence of the event rates in COHERENT, and the statistical analysis used. Our results are then presented in section 4, first using only COHERENT data (section 4.1) and then in combination with the results from oscillation experiments (section 4.2). In doing so we discuss the dependence of the results on variations of the COHERENT analysis associated to the choice of quenching factor, the nuclear form factor, and the treatment of the backgrounds. Finally in section 5 we summarize and draw our conclusions.

Coherent elastic neutrino-nucleus scattering
At the SNS, an abundant flux of both π + and π − is produced in proton-nucleus collisions in a mercury target. While the π − are absorbed by nuclei before they can decay, the π + lose energy as they propagate and eventually decay at rest into π + → µ + ν µ , followed by µ + → e + ν eνµ . Since the muon lifetime is much longer than that of the pion, the ν µ component is usually referred to as the prompt contribution to the flux, as opposed to the delayed contributions from µ + decay (ν µ and ν e ).
Given that the prompt neutrinos are a by-product of two-body decays at rest, their contribution to the total flux is a monochromatic line at E pr = (m 2 π − m 2 µ )/(2m π ) 29.7 MeV, where m π and m µ refer to the pion and muon masses, respectively. Conversely, the delayed neutrino fluxes follow a continuous spectra at energies E νe,νµ < m µ /2 52.8 MeV. At a distance from the source, they read: and are normalized to each proton collision on the target. For reference the distance at COHERENT is 19.3 m. The differential cross section for coherent elastic neutrino-nucleus scattering, for a neutrino with incident energy E ν interacting with a nucleus with Z protons and N neutrons, reads [34]:

JHEP02(2020)023
where G F is the Fermi constant and Q 2 is the weak charge of the nucleus. In this notation, T is the recoil energy of the nucleus, M is its mass, and F is its nuclear form factor (FF) evaluated at the squared momentum transfer of the process, Q 2 = 2M T . In our calculations we have first used a Helm FF 2 parametrization [51]: where s = 0.9 fm [52] and j 1 (x) is the order-1 spherical Bessel function of the first kind. The value of R 0 relates to the value of s and the neutron radius R n as In the absence of an experimental measurement of the neutron radius in CsI, we tune its value so that the prediction for the total number of events at COHERENT matches the official one provided in refs. [35,36] (173 events). However, by doing so we obtain R n = 4.83 fm, a value that is unphysical as it approaches the proton radius [53] which all models predict to be smaller. Given that this is a phenomenological parametrization, though, and seeing the large differences in the prediction for the total number of events obtained with different values of R n , it is worth asking whether this is accurate enough for CsI, and exploring the impact of the nuclear FF on the results of the fit. Therefore, in section 4 we will also show the results obtained using a state-of-the-art theoretical calculation for the nuclear FF, taken from refs. [48,49] (calculated using the same methodology as in refs. [54,55]). In the SM, the weak charge of a nucleus only depends on the SM vector couplings to protons (g V p ) and neutrons (g V n ) and is independent of the neutrino flavor: where g V p = 1/2 − 2 sin 2 θ w and g V n = −1/2, with θ w being the weak mixing angle. For CsI, we obtain Q 2 1352.5 in the SM. However, in presence of NC NSI, this effective charge gets modified 3 by the new operators introduced as [56]: and in general its value may now depend on the NSI parameters ε ≡ {ε f αβ } as well as the incident neutrino flavor α. Since the COHERENT experiment observes interactions of both electron and muon neutrinos, its results are sensitive to both Q 2 e and Q 2 µ .

JHEP02(2020)023
As can be seen from eq. (2.6), the modification of NSI to the CEνNS event rate comes in as a normalization effect. Therefore, adding energy information to the analysis of the data is not expected to have a significant effect on the results of our fit to NSI. 4 Conversely, the addition of timing information is crucial as it translates into a partial discrimination between neutrino flavors, thanks to the distinct composition of the prompt (ν µ ) and delayed (ν µ and ν e ) neutrino flux. This translates into an enhanced sensitivity to NSI, since the fit will now be sensitive to a change in normalization affecting neutrino flavors differently.

Implementation of the COHERENT experiment
In our previous work [29] we performed a fit to COHERENT data using the available information at that time, which included only the total event rates observed. Last year the collaboration released publicly both the energy and timing information of the events [36]. In this section we describe the procedure used to implement in our fit the information provided in such data release.

Computation of the signal
The differential event distribution at COHERENT, as a function of the nuclear recoil energy T , reads dN dT where N nuclei is the total number of nuclei in the detector, N pot = 1.76 · 10 23 is the total number of protons on target considered, and f ν/p = 0.08 is the neutrino yield per proton. In eq. (3.1) the sum runs over all neutrino flux components (ν e , ν µ ,ν µ ), and the upper limit of the integral is given by the end-point of the spectrum from pion DAR, while the minimum neutrino energy that can lead to an event with a nuclear recoil energy T is given by At COHERENT, the observable that is actually measured is the number of photoelectrons (PE) produced by an event with a certain nuclear recoil. In fact, the nuclear recoil energy in CEνNS events is typically dissipated through a combination of scintillation (that is, ionization) and secondary nuclear recoils (that is, heat). While secondary recoils are the characteristic signal of a nuclear recoil (as opposed to an electron recoil, which favors ionization instead), their measurable signal is much smaller than that of electron recoils. The ratio between the light yields from a nuclear and an electron recoil of the same energy is referred to as the Quenching Factor (QF).
Besides being a detector-dependent property, the QF may also depend non-trivially on the recoil energy of the nucleus. In general, the relation between PE and nuclear recoil T

JHEP02(2020)023
can be expressed as: where LY is the light yield of the detector (that is, the number of PE produced by an electron recoil of one keV), and we have explicitly noted that the QF may depend on the nuclear recoil energy. Therefore, the expected number of events in a certain bin i in PE space can be computed as: where the limits of the integral correspond to the values of T obtained for the edges of the PE bin (PE min i , PE max i ) from eq. (3.3). In their analysis, the COHERENT collaboration adopted an energy-independent QF throughout the whole energy range considered in the analysis, between 5 and 30 keV [35]. Also, given the tension observed between the different calibration measurements available at the time, they assigned large error bars to the assumed central value QF = 8.78%. Taking a central value for the light yield LY = 13.348 PE per keV of electron recoil [36], this means that approximately 1.17 PE are expected per keV of nuclear recoil energy. Very recently, however, the authors of ref. [58] have re-analyzed past calibration data used to derive this result. They concluded that the tension between previous measurements was partially due to an unexpected saturation of the photo-multipliers used in the calibration and, after correcting for this effect, a much better agreement was found between the different data sets. This allows for a significant reduction of the error bars associated to the QF, as well as for the implementation of an energy-dependent QF. On the other hand, the COHERENT collaboration has not confirmed the claims of the authors of ref. [58]. After repeating the calibration measurements, their new data still shows a good agreement [59] with the original measurements performed by the Duke (TUNL) group [35,36].
Given that this issue has not been settled yet, we will study and quantify the effect of these new measurements in the results of the fit in section 4. We will present our results obtained for three different QF parametrizations: the original (constant) parametrization used in the data release [36]; the best-fit obtained by the authors of ref. [58]; and the results from our fit to the calibration measurements performed by the TUNL group [35,36]. In order to fit the data of the Duke group, we use the phenomenological parametrization proposed in [58], which is based on a modification of the semi-empirical approach by Birks [60] and depends only on two parameters, E 0 and kB. 5 Following this approach, we obtain a best-fit to the Duke group data for E 0 = 9.54 ± 0.84 and kB = 3.32 ± 0.10, with a correlation ρ kB,E 0 = −0.69. The three QF parametrizations used in our calculations are shown in figure 1, as a function of the recoil energy of the nucleus.
Once the expected event distribution in PE has been computed following eq. , provided as part of the data release [36]. In both panels, the constant QF used in ref.
[35] is also shown for comparison. The three parametrizations are shown with a shaded band to indicate the values allowed at the 1σ CL in each case. For illustration, the vertical shaded area indicates approximately the range of nuclear recoil energies that enters the signal region used in the fit (the exact range varies with the nuisance parameters, and the exact QF parametrization used). distribution, to account for the probability that a given event yields a different number of PE than the average. On top of that, signal acceptance efficiencies are applied to each bin: where the function Θ is defined as: Following ref. [36], the central values of the signal acceptance parameters are set toη 0 = 0.6655,k = 0.4942, PE 0 = 10.8507. Finally, once the predicted energy spectrum has been computed, one should consider the arrival times expected for the different contributions to the signal. This is implemented using the distributions provided by the COHERENT collaboration in the data release [36], which are normalized to one.
This final prediction can be compared with the published data. This is provided in two different time windows for each trigger in the data acquisition system (that is, for each proton pulse). On the one hand, the region where signals and beam-induced backgrounds associated with the SNS beam are expected is referred to as the coincidence (C) region, which can therefore be considered a "signal" region. Conversely, the region where no contribution from the SNS beam is expected is referred to as the anti-coincidence (AC) region and could be considered a "background" region. While the collaboration provides data separately for the beam-ON and beam-OFF data taking periods, in this work we JHEP02(2020)023 only use the beam-ON samples. The total exposure considered in this work corresponds to 308.1 live-days of neutrino production, which correspond to 7.48 GW-hr (∼ 1.76 × 10 23 protons on target). The residual event counts for this period, i.e., the C data with the AC data subtracted, are shown in figure 2, projected onto the time and PE axes, for different choices of the QF and FF as indicated by the labels.
As can be seen from the comparison between the upper and lower panels, the change in QF between a constant approximation (Data Release) and the energy-dependent result obtained by the TUNL group (QF-D) does not affect significantly the predicted event distributions. This will lead to a minor change in the results of the numerical fit to the data in section 4. A larger difference is observed with respect to the predictions using the QF by the Chicago group (QF-C, middle panels): in this case, the very different central values at T ∼ 10 keV (corresponding to PE ∼ 10) lead to a reduced number of events, which will have a larger impact on the results.

Computation of the background
The COHERENT measurement is affected by three main background sources: (i) the steady-state background, coming from either cosmic rays or their by-products entering the detector; (ii) prompt neutrons produced in the target station and exiting it, and (iii) neutrino-induced neutrons (NINs) that originate in the shielding surrounding the detector. While the latter is irreducible, it has been shown to be negligible at the COHERENT experiment and is therefore ignored here.
The procedure used to compute the expected number of background events for the steady-state and the prompt neutron components follows the prescription given in ref. [36]. For both backgrounds, it is assumed that the temporal and energy dependence on the number of events can be factorized as where f contains the temporal dependence of the signal and g its energy dependence. For the prompt neutron background, the collaboration provides both its expected energy distribution before acceptance efficiencies are applied (that is, g(PE)), and the total expected counts as a function of time (that is, f (t)). The expected 2D distribution can be obtained simply by multiplying the two distributions. After the number of events in each bin has been computed, the same acceptance efficiency as for the signal, eq. (3.5), is applied to determine the expected number of events in each bin.
The steady-state is the most significant background source to this analysis, and it has the largest impact on the fit. In this case, the functions g(PE) and f (t) are not provided in ref. [36] but inferred from the data, which is provided per bin in energy and time. In particular, the projected data onto the PE axis is then used directly as g(PE), while f (t) is assumed to follow an exponential: In the upper panels the shaded histograms show the predicted event rates in the SM using the QF and nuclear FF from ref. [36]. In the middle panels they correspond to the predictions with the QF from the Chicago group (QF-C) in ref. [58] and the nuclear FF from ref. [48,49]. The lower panels have been obtained with the same FF, but changing the QF to match the Duke (TUNL) measurements in ref.
[35] (QF-D). In all panels the prompt neutron background prediction is also shown for completeness. All the event histograms shown in this figure correspond to the SM prediction. so that its integral over the whole range in time is equal to one. Since in this case the expected background events are inferred from a measurement, the signal acceptance has already been included into the calculation and there is no need to apply it here.
JHEP02(2020)023 The procedure outlined above for the steady-state component is meant to eliminate biases in the fit due to the limited statistics of the data sample used. However, by treating the background in this way the analysis is rather sensitive to a mis-modeling of its temporal component. In particular, if we plot the separate C and AC event distribution as a function of time, instead of looking at their difference, it is easy to see that there is an excess in the first two bins in the data, which cannot be accommodated by the simplified exponential fit. This is shown in figure 3, where we show the total AC counts (which should include only the steady-state background as measured by the detector) together with the exponential that gives a best fit to the data. As clearly seen, the first two bins are not well fitted by a simple exponential model.
Interestingly enough, both the C and AC samples seem to observe a similar excess in the first two temporal bins, which suggests that this contribution is not related to the neutrino signal but to some mis-modeling of the background. A possible way to correct for this is to directly use the measured time dependence of the AC sample as a direct prediction for the expected behavior of the steady-state background in the C sample. Doing this on a bin-per-bin basis would not provide a good predictor for the expected number of events in each bin, due to the limited statistics. However, the projected data onto the time axis may still be used, as in the case of the exponential fit, to get a prediction for the function f (t). In other words, the prediction for f (t) may be obtained following the same procedure as was done for g(PE), i.e., projecting the events onto the corresponding axis. Figure 4 shows the observed total event counts for the beam-ON, C sample (which includes both signal and background) projected onto PE (left) and time (right), compared with the predictions using these two different background models.
As can be seen from this figure, both background models are able to reproduce the observed spectrum in PE relatively well, and give very similar results. However, the event rates obtained with this second method provide a better fit to the data when projected JHEP02(2020)023   [36], while the lower panels have been obtained using our model (which follows the time dependence of the AC events, see text for details). The light histograms show the predicted total number of events, after adding all signal and background contributions. To ease the comparison among different panels, in this figure the signal has been computed in all cases using the same FF and QF as in the data release [36]. All histograms shown correspond to the SM predicted event rates.
onto the time axis and, as clearly observed from the figure, the effect is specially noticeable in the first two bins. Therefore, in section 4 we will show two sets of results: with and without using an exponential model for the background.

Systematic errors and implementation of the χ 2
Once the predicted event distributions for the signal and backgrounds have been computed, a χ 2 function is built as: where O ij stands for the observed number of events in PE bin i and time bin j, while P ij stands for the total number of predicted events in that bin, including the signal plus all background contributions. Following ref. [36], we consider only the events with 5 < PE ≤ 30

JHEP02(2020)023
and t < 6 µs in the analysis. The predicted number of events depends on the nuisance parameters ξ ≡ {ξ a } included in the fit, which account for the systematic uncertainties affecting the QF, signal acceptance, neutrino production yield, and normalization of the backgrounds. These are implemented replacing the original quantity as x → (1 + σ x ξ x )x, wherex denotes the central value assumed for x prior to the experiment and σ x denotes the relative uncertainty for nuisance parameter ξ x summarized in table 1 for convenience. More specifically: where N ss ij , N n ij and N sig ij stand for the predicted number of events for the steady-state background, the prompt neutron background and the signal. In eq. (3.10) we have generically denoted as ξ QF the set of nuisance parameters characterizing the uncertainty on the QF employed. For the constant parametrization used in the data release [36] we introduce a unique nuisance parameter with constant uncertainty. For QF-C we introduce also a unique nuisance parameter, but with an energy-dependent uncertainty inferred from the uncertainty band in figure 1 of ref. [58] (also shown in the left panel of figure 1), which varies from 6.5% to 3.5% in the range of recoil energies relevant for COHERENT. For QF-D we introduce two nuisance parameters characterizing the uncertainty on parameters E 0 and kB with their corresponding correlation. 6 Altogether the likelihood for some physics model parameters ε, leading to a given set of predictions P ε ij ( ξ) for the events in bin ij, is obtained including the effects of the nuisance parameters as in eq. (3.10) and minimizing over those within their assumed uncertainty. This is ensured by adding a pull term to the χ 2 function in eq. (3.9) for each of the nuisance parameters introduced: where ρ is the correlation matrix, whose entries are ρ ab = δ ab for all parameters except those entering our parametrization of the QF of the Duke group (the corresponding correlation coefficient can be found in table 1). As validation of our χ 2 construction we have performed a fit to extract the total number of signal CEνNS events when using the same assumptions on the background, systematics and energy and time dependence of the signal as those employed by COHERENT in their data release [36]. This can be directly compared with their corresponding likelihood extracted from figure S13 of ref. [35]. The result of this comparison is shown in figure 5. Strictly speaking, the χ 2 function plotted in figure 5 depends on the assumed energy dependence of the signal. Therefore it is expected to vary if, instead of using the QF and FF quoted in the data release, we employed a different QF parametrization and nuclear FF.
5.0 Prompt n norm. 25.0 Signal norm. 11.2 η 0 4.5 k 4.7 PE 0 2.7 QF (data release) 18.9 QF (our fit, QF-C) 6.5-3.5 kB (our fit, QF-D) 3.0 E 0 (our fit, QF-D) 8.8 ρ E 0 ,kB (our fit, QF-D) −0.69 Table 1. Systematic uncertainties considered in the fit on acceptance efficiency parameters (eq. (3.5)), normalization of the signal and background contributions, and the QF. The steadystate normalization uncertainty includes the statistical error of the sample (AC data). The quoted uncertainties on the QF also includes the error on the light yield (0.14%), which is however subdominant. For details on the QF parametrization, see section 3.
Quantitatively, within the systematic uncertainties used in the construction of the shown χ 2 (N CEνNS ), we find that changing the QF and FF has a negligible effect on this curve. Conversely, we find a stronger dependence on the systematic uncertainties introduced, and therefore figure 5 serves as validation of our implementation for these. For illustration, we also indicate the predicted event rates in the SM predicted by the collaboration (173 events) as well as our result obtained using the Chicago QF parametrization [58], as shown in the left panel of figure 1 and the new nuclear FF from refs. [48,49]. In both cases the vertical lines correspond to the prediction without accounting for systematic uncertainties. The predicted result using the Duke QF and the new FF from refs. [48,49] is very similar to the one obtained by the collaboration (168 events) and is therefore not shown here.

Results
In this section we present our results. First, in section 4.1 we provide the results of the new fit using COHERENT data alone, which we compare with those obtained in our previous work [30]. We discuss in detail the improvements coming from the inclusion of energy and timing information in the fit, and investigate the impact of the different choices of quenching factor, nuclear form factor and background implementation. We then proceed to combine the COHERENT data with the global analysis of oscillation data in section 4.2. We will show two main sets of results, which quantify: (1) the quality of the global fit once NSI are allowed, compared to the fit obtained under the SM hypothesis; and (2) the status of the LMA-D solution after the inclusion of COHERENT energy and timing data in the global fit. Both sets of results are presented for a wide range of NSI models (that is, for different values of η). Finally, we also provide the allowed JHEP02(2020)023  [36], while solid green indicates our prediction using the QF from ref. [58] (left panel in figure 1) and nuclear FF of [48,49]. ranges obtained for the Wilson coefficients for three particular NSI models, assuming that the new mediator couples predominantly to either up/down quarks or to protons.

Fit to COHERENT data
In order to study the dependence of the results on the different assumptions, we have performed a set of fits to COHERENT data in terms of two effective flavor-dependent weak charges Q 2 α (assumed to be energy-independent). Figure 6 shows the corresponding allowed regions at 2σ from the fit to COHERENT data alone. In the left panel we illustrate the effect of including the energy and timing information in the fit, by comparing the allowed values of the weak charges obtained: (i) using only the total rate information (dot-dashed); (ii) adding only the energy information (dashed); (iii) using only the event timing information (dotted); and (iv) fitting the data binned in both timing and energy (solid). In all cases shown in the left panel, the QF, nuclear FF and background have been implemented following closely the prescription given in the data release [36]. It is wellknown that, when only the total event rate information is considered, there is a degeneracy in the determination of the flavor-dependent weak charges, since the number of predicted events approximately behaves as where f α indicates the fraction of expected SM events from interactions of ν α in the final event sample and we have assumed that one neutrino of each species is produced for each pion DAR. Under these assumptions, the allowed region in the (Q 2 e , Q 2 µ ) plane is a straight band with a negative slope, arctan(−0.5) ≈ −27 • . This behavior is also observed from our exact fit to the data, as shown by the dot-dashed lines in the left panel in figure 6.
As expected, the timing information is most relevant in breaking of this degeneracy: since the prompt component of the beam contains only ν µ , the inclusion of time information allows for a partial discrimination between Q 2 e and Q 2 µ . Notice, however, that in this case the best fit is obtained at the edge of the physically allowed region, Q 2 e 0 (in fact, it would probably take place for a negative value, but this is not the case since we are effectively imposing the restriction Q 2 α > 0 in the fit). This is driven by the small excess for the event rates in the first two time bins (with respect to the SM prediction) when using the exponential fit model for the steady-state background, as described in section 3.2 (see figure 4). Such excess can be accommodated thanks to the overall normalization uncertainty of the signal, combined with a decrease of the ν e contribution as required to match the distribution observed for the delayed events. Within the systematic uncertainties in the analysis, this results into a higher rate at short times without a major distortion of the PE spectrum. We also observe that, including only the PE spectrum in the fit, the degeneracy still remains but the width of the band in this plane decreases. For values of Q 2 α in the non-overlapping region, the fit using the event rate information alone is able to fit the data, albeit at the price of very large nuisance parameters and, in particular, of the QF-related uncertainties (which affect the shape of the event distributions in PE space). Therefore, once the PE information is added the allowed regions are consequently reduced.  Figure 7. 2σ allowed regions for the flavor-diagonal NSI coefficients ε coh αβ (assuming zero nondiagonal couplings) for a variety of fits to COHERENT data as labeled in the figure. In all cases shown in the left panel, the QF, nuclear FF and background assumptions are those employed in the data release [36]. On the right panel we show the dependence on the assumptions for steady background modeling, nuclear FF, and QF. For simplicity, in this figure we set all off-diagonal NSI parameters to zero, but it should be kept in mind that the results of our global analysis presented in section 4.2 have been obtained allowing all operators simultaneously in the fit.
The right panel in figure 6 shows the dependence of the allowed region on the assumed background model, nuclear FF and QF choice in the fit, for the 2D fit using both time and PE information. As seen in the figure, if one uses the steady-state background prediction without the exponential model for its temporal dependence the region becomes considerably larger, and the BF moves closer to the SM. This is expected because, with this background, there is no excess of events in the first time bins with respect to the SM prediction (see figure 4). This also leads to a better overall fit, with χ 2 min = 145.24 (for 12 × 12 = 144 data points) compared to χ 2 min = 150.8 obtained for the exponential model of the steadystate background. Altogether we observe that modifying the nuclear FF has a very small effect on the current results, as can be seen from the comparison between the orange and green lines in the figure. Changing the QF does not have a dominant impact either, once the exponential fit to the background has been removed. This can be seen from the comparison between the green, brown and blue lines in the figure, which all provide similar results. Overall, we find a slightly better agreement with the SM result for the QF-C parametrization, albeit the effect is small.
In the framework of NSI, the constraints on the weak charges derived above can be directly translated into constraints on the effective Wilson coefficients. This is shown in figure 7, where we plot the allowed regions for the two relevant flavor-diagonal NSI couplings after setting the flavor-changing ones to zero. In doing so we notice that, as JHEP02(2020)023 discussed in ref. [28], the fact that the neutron/proton ratio in the two target nuclei is very similar (N Cs Z Cs 1.419 for cesium and N I Z I 1.396 for iodine) allows to approximate eq. (2.6) as: with an average value Y coh n = 1.407 and From eq. (2.6) it is evident that COHERENT can only be sensitive to a certain combination of NSI operators ε coh αβ , which are ultimately determined by just two factors: (a) the value of Y coh n , which depends on the nuclei in the detector, and (b) the strength of the coupling of the new interaction to up and down quarks (or, equivalently, to protons and neutrons). In fact, using the η parametrization in eq. (1.4), ε coh αβ can be written as: It is clear from the expressions above that the best-fit value and allowed ranges of ε coh αβ implied by COHERENT are independent of η. Once these have been determined, the corresponding bounds on the associated couplings ε η αβ for a given NSI model (identified by a particular value of η) can be obtained in a very simple way, by just rescaling the values of ε coh αβ as [ √ 5(cos η + Y coh n sin η)] −1 . For example, the results in figure 7 can be immediately translated in the corresponding ranges for NSI models where the new interaction couples only to f = u, f = d, or f = p (η ≈ 26.6 • , 63.4 • , and 0, respectively), after rescaling the bounds on ε coh αβ by the corresponding factors of 0.293, 0.262, and 1 in each case. Furthermore, from eq. (4.4) it becomes evident that, for NSI models with η = arctan(−1/Y coh n ) ≈ −35.4 • , no bound can be derived from COHERENT data. The impact of the timing information on the fit can be readily observed from the left panel in figure 7. In the absence of any timing information and using total rate information alone, it is straightforward to show that, if the experiment observes a result compatible with the SM expectation, the allowed confidence regions in this plane should obey the equation of an ellipse. This automatically follows from eqs. (4.1) and (4.2): where R ≡ g V p + Y coh n g V n ≈ −0.68. This is also shown by our numerical results in the left panel of figure 7 which do not include timing information in the fit (dashed and dot-dashed contours).
While the inclusion of a non-zero ε coh ee can be compensated by a change in ε coh µµ that brings the total number of events in the opposite direction without significantly affecting the delayed events, this would be noticed in the prompt event distribution once timing information is added to the fit. In particular, too large/small values of ε coh ee would require a consequent modification of the ε coh µµ to recover the same event rate in the delayed time bins, which is however not allowed by the prompt events observed. Thus, once timing JHEP02(2020)023 information is included in the fit the ellipse is broken in this plane and two separate minima are obtained (dotted and solid lines).
It should also be noted that the central region in figure 7 (around the centre of the ellipse in eq. (4.5), ε coh ee = ε coh µµ = −R) can be excluded at COHERENT only in the case when the off-diagonal NSI operators are not included in the fit. This is so because in this region the effect of the diagonal parameters leads to a destructive interference in the total cross section and therefore to a reduction of the number of events, in contrast with the experimental observation. Once the off-diagonal operators are introduced this is no longer the case and the central region becomes allowed [39]. However, since global neutrino oscillation data provide tight constraints on the off-diagonal NSI operators, in our results the two minima remain separate even after the off-diagonal operators are allowed in the fit, as we will show in section 4.2.

Global analysis of COHERENT and oscillation data
We now present the results of the global analysis of oscillation plus COHERENT data. To this end we construct a combined χ 2 function where we denote by ω ≡ {θ ij , δ CP , ∆m 2 ji } the "standard" 3ν oscillation parameters. For the detailed description of methodology and data included in χ 2 OSC we refer to the comprehensive global fit in ref. [28] performed in the framework of three-flavor oscillations plus NSI with quarks parametrized as eq. (1.4). In this work we minimally update the results from ref. [28] to account for the latest LBL data samples included in NuFIT-4.1 [62]. To keep the fit manageable in ref. [28] only the CP-conserving case with real NSI and δ CP ∈ {0, π} was considered, and consequently the T2K and NOνA appearance data (which exhibit substantial dependence on the leptonic CP phase) were not included in the fit. Here we follow the same approach and consistently update only the disappearance samples from these experiments. 7 Figure 8 shows the impact of COHERENT on the global fit (left panel) as well as on the LMA-D degeneracy (right panel). In doing so, we have defined the functions χ 2 LMA (η) and χ 2 LMA-D (η), obtained by marginalizing χ 2 global ( ω, ε) over both ω and ε for a given value of η, with the constraint θ 12 < 45 • (in the LMA case) and θ 12 > 45 • (for the LMA-D). With these definitions, we show in the left panel the differences χ 2 LMA (η) − χ 2 no-NSI (full lines) and χ 2 LMA-D (η) − χ 2 no-NSI (dashed lines), where χ 2 no-NSI is the minimum χ 2 for standard 3ν oscillations (i.e., setting all the NSI parameters to zero). Then, in the right panel we show the values of χ 2 LMA-D (η) − χ 2 LMA (η), which quantifies the relative quality of the LMA and LMA-D solutions as a function of η.
First, from the left panel in figure 8 we notice that the introduction of NSI leads to a substantial improvement of the fit already for the LMA solution (solid lines) with respect to the oscillation data analysis, resulting in a sizable decrease of the minimum χ 2 LMA with respect to the standard oscillation scenario. This is mainly driven by a well-known tension 7 For a discussion of CP violation in the presence of NSI see ref. [63]. (although mild, at the level of ∆χ 2 ∼ 7.4 in the present analysis) between solar and KamLAND data in the determination of ∆m 2 21 . As seen in the figure, the inclusion of NSI improves the combined fit by about 2.2σ over a broad range of values of η. The improvement is maximized for NSI models with values of η for which the effect is largest in the Sun without entering in conflict with terrestrial experiments. This occurs for η −44 • (as for this value the NSI in the Earth matter essentially cancel) and leads to an improvement of about 10 units in χ 2 (i.e., a ∼ 3.2σ effect). From the figure we also conclude that adding the information from COHERENT on rate only, as well as on timing and energy (t+E), still allows for this improved fit in the LMA solution for most values of η. Indeed, the maximum effect at η −44 • still holds after the combination since it falls very close to −35.4 • , for which NSI effects cancel at COHERENT as seen in eq. (4.4). Interestingly, the improvement is slightly larger for the combination with COHERENT t+E data using the data release assumptions. This is so because, as described in the previous section, in this case the fit pulls the weak charge Q 2 e towards zero (see figure 6) while leaving the value of Q 2 µ around the SM expectation. Such situation can be easily accommodated by invoking diagonal NSI operators and, in particular, favors the non-standard values ε η ee − ε η µµ = 0, thus bringing the fit to a better agreement with solar+KamLAND oscillation data.
Most importantly, figure 8 shows that the main impact of including COHERENT data in the analysis is on the status of the LMA-D degeneracy. We see in the figure that with JHEP02(2020)023 oscillation data alone the LMA-D solution is still allowed at 3σ for a wide range of NSI models (−38 • η 87 • , as well as a narrow window around η −65 • ) and, in fact, for −31 • η 0 • it provides a slightly better global fit than the LMA solution. The addition of COHERENT to the analysis of oscillation data disfavors the LMA-D degeneracy for most values of η, and the inclusion of the timing and energy information makes this conclusion more robust. More quantitatively we find that, when COHERENT results are taken into account, the LMA-D is allowed below 3σ only for values of η in the following ranges: Finally, we provide in figure 9 the χ 2 profiles for each of the six NSI coefficients after marginalization over the undisplayed oscillation parameters and the other five NSI coefficients not shown in a given panel. We show these results for three representative cases of NSI models including couplings to up quarks only, down quarks only and to protons. The corresponding 2σ ranges are also provided in table 2 for convenience. This figure shows that the LMA-D solution for NSI models that couple only to protons (η = 0) can only be excluded beyond 3σ if both energy and timing information are included for COHERENT, in agreement with eq. (4.7). From this figure we also see that for the LMA solution the allowed ranges for the off-diagonal NSI couplings are only moderately reduced by the addition of the COHERENT results and, moreover, the impact of the energy and timing information is small in these cases. This is because they are already very well constrained by oscillation data alone.
More interestingly, the addition of COHERENT data allows to derive constraints on each of the diagonal parameters separately and, for those, the timing (and to less degree energy) information has a quantitative impact. In particular we see in the figure the appearance of the two minima corresponding to the degenerate solutions for ε coh µµ in figure 7, obtained after the inclusion of timing information for COHERENT. Noticeably, now the non-standard solution (obtained for ε f µµ = 0) is partially lifted by the combination with oscillation data, but it remains well allowed around ∼ 2σ depending on the assumptions for the COHERENT analysis. Figure 9 also shows the corresponding two minima for ε f τ τ arising from the combination of the information on ε f τ τ − ε f µµ and ε f ee − ε f µµ from the oscillation experiments with the constraints on ε f µµ from the COHERENT t+E data. In particular, in all the three cases f = u, d, p shown in the figure the bound on ε f τ τ becomes about two orders of magnitude stronger than previous indirect (loop-induced) limits [5] when the t+E analysis with the data release assumptions is used. Indeed this conclusion holds for most η values, with exception of η ∼ −45 • to −35 • for which NSI effects are suppressed in either the Earth matter or in COHERENT. In particular, for η = −35.4 • the NSI effects in COHERENT totally cancel, as described above, and consequently no separate determination of the three diagonal parameters is possible around such value.  . Dependence of the ∆χ 2 global function on the NSI couplings with up quarks (upper row), down quark (central row) and, protons (lower row) for the global analysis of oscillation and COHERENT data. In each panel χ 2 global is marginalized with respect to the other five NSI couplings not shown and with respect to the oscillation parameters for the LMA (solid) and LMA-D (dashed) solutions. The different curves correspond to the different variants of the COHERENT analysis implemented in this work: total rate (black), t+E Data Release (red), t+E with QF-C (blue), and t+E with QF-D (brown); see text for details.

Summary and conclusions
From a completely model-independent approach, a useful way to parametrize the effects of new physics models on low-energy experiments is through the inclusion of additional higher-dimensional operators to the Standard Model Lagrangian. In this context, fourfermion operators leading to neutral-current interactions between neutrinos and quarks, usually referred to as neutrino non-standard interactions (NSI), can lead to novel effects in both neutrino propagation, production and detection processes. Although models leading to large NSI in the neutrino sector are challenging to build from a high-energy theory, JHEP02(2020)023 new physics models invoking light mediators are able to evade the tight constraints in the charged lepton sector, and can lead to observable consequences in neutrino experiments. While charged-current NSI are better constrained by precision measurements of muon and meson decays, neutral-current NSI are harder to constrain experimentally and, in fact, the best bounds in this case come from global fits to neutrino data.
In this work, we have combined neutrino oscillation data with the latest results obtained for coherent neutrino-nucleus scattering data at the COHERENT experiment, which provide both energy and timing information. The results of our analysis are used to constrain the Wilson coefficients of the whole set of neutral-current operators leading to NSI involving up and down quarks simultaneously. Therefore, our conclusions extend to a general class of new physics models with arbitrary ratios (parametrized by an angle η) between the strength of the operators for up and down quarks, and mediated by Z with masses above M med ∼ O(10) MeV.
We have quantified the dependence of our results for COHERENT with respect to the choice of quenching factor, nuclear form factor, and the treatment of the backgrounds. We find that the implementation of the steady-state background has a strong impact on the results of the analysis of COHERENT due to a slight background excess in the first two bins, which is present in both the coincident and anti-coincident data samples provided by the collaboration. Once this effect has been accounted for in the modeling of the expected backgrounds, the choice of quenching factor and nuclear form factor has a minor impact on the results obtained from the fit.
We find that the inclusion of COHERENT timing information affects the global fit significantly and, most notably, has a large impact on the constraints that can be derived for the flavor-diagonal NSI operators, for which separate constraints can only be derived after combination of COHERENT and oscillation data.
Furthermore, the presence of NSI is known to introduce a degeneracy in the oscillation probabilities for neutrinos propagating in matter, leading in particular to the appearance of the LMA-Dark (LMA-D) solution. We find that the inclusion of COHERENT to the analysis of oscillation data disfavors the LMA-D degeneracy for NSI models over a wide range of η, and the addition of the timing and energy information makes this conclusion more robust, see eq. (4.7) and figure 8. In particular, the LMA-D solution for NSI models that couple only to protons (η = 0) can only be excluded beyond 3σ once both energy and timing information are included for the COHERENT data.
Finally, the introduction of NSI is known to alleviate the well-known (albeit mild) tension between solar and KamLAND data in the determination of the ∆m 2 21 , thus leading to an overall improvement in the quality of the global fit to oscillation data. We find that this still remains the case after the inclusion of COHERENT results.