Searching for Beyond the Standard Model Physics with COHERENT Energy and Timing Data

We study the prospects for extracting information on beyond the Standard Model (BSM) physics by combining the energy and timing distribution of nuclear recoil events in the COHERENT neutrino-nucleus elastic scattering data. Focusing on light, $\lesssim$ GeV mediators as an example, we find that the combined energy and timing data favor a $\sim 10-1000$ MeV mediator as compared to the Standard Model (SM) best fit at the $\lesssim 2\sigma$ level. The timing data provides additional statistical information on the electron and muon neutrino flavor distributions that is not attainable from the energy distribution of nuclear recoils alone. This result accounts for the uncertainty in the effective size of the neutron distribution, and highlights the power of an analysis which includes uncertainties on the experimental backgrounds, the nuclear model, the BSM parameters, and uses the full energy and timing information.


I. INTRODUCTION
The COHERENT collaboration has reported the first detection of coherent neutrinonucleus elastic scattering (CEνNS) [1].COHERENT utilizes the Spallation Neutrino Source (SNS) with a stopped-pion beam, which produces a well-known neutrino spectrum from pion and muon decay at rest.Muon neutrinos, ν µ , thus arrive from a prompt decay of charged pions, whereas νµ and ν e are produced from the delayed decay of muons.With an exposure of 14.6 kg-days, the COHERENT collaboration measured a best fit of 134 ± 22 nuclear recoil events from CEνNS, which is well in excess of the expected background events for this exposure.
The COHERENT data provides an important new channel to search for beyond the Standard Model (BSM) physics.For example, the COHERENT data constrains non-standard neutrino interactions (NSI) [2,3] which may arise from heavy or light mediators [4][5][6][7][8].The data also provide novel constraints on new physics that manifests through the neutrino sector, including generalized scalar and vector neutrino interactions [9], as well as hidden sector models [10].It also sets independent constraints on the effective neutron size distribution of CsI [11][12][13], and on sterile neutrinos [14,15].Low mass mediator models are particularly interesting since there are no Large Hadron Collider (LHC) constraints for mediator masses < ∼ GeV, and these models connect to new ideas associated with sub-GeV dark matter in the range MeV to GeV [16].The aforementioned constraints on BSM physics have mostly been extracted from the energy distribution of nuclear recoil events.
The COHERENT data do not directly identify the flavor components of the neutrino flux, though it is possible to make an estimate for the contribution of the different flavors.
For example the prompt ν µ component may be estimated from a timing cut, while the delayed components may be extracted from their spectral signatures.Previous authors have classified the events as prompt or delayed using a two-bin analysis in timing space [17].
Including timing information in this manner has the potential to strengthen constraints on NSI parameters, in comparison to those obtained from using information in the energy distribution alone.
In this paper, we perform a likelihood analysis on the COHERENT data that utilizes the full energy and timing distributions of nuclear recoil events.These distributions provide information on the flavor components of the detected neutrino flux that is not available when considering the energy data alone, or when splitting the timing data into prompt and delayed events.We test for deviations from pure Standard Model (SM) interactions, considering as an example light mediators that couple to the SM.We also consider a scenario in which the quark-Z coupling contains a momentum dependent factor associated with hidden sector interactions.We show that including the timing distribution data adds substantial statistical constraining power, and in fact favors a BSM interpretation of the nuclear recoil data at the < ∼ 2σ level.

II. LIGHT MEDIATORS AND NON-STANDARD INTERACTIONS
In this section we briefly discuss the SM prediction for CEνNS.We introduce a model for light mediators as an example of BSM physics, and highlight how this model alters the predicted CEνNS cross section.
In the SM, the differential cross section for a neutrino scattering off of a target electron or quark of mass m through a Z exchange is where G F is the Fermi constant, E is the recoil energy, E ν is the incident neutrino energy, (g v , g a ) = (T 3 − 2Q em sin 2 θ W , T 3 ) are the vector and axial-vector couplings to the Z-boson, T 3 is the third component of the weak isospin, Q em is the electromagnetic charge, and θ W is the weak mixing angle (T 3 = −1/2 in our convention).For quarks in the nucleus, if the momentum transfer between the neutrino and the nucleus is smaller than or comparable to the inverse size of the nucleus (typically E < ∼ O(50) MeV), coherent scattering occurs.For coherent elastic scattering off of nuclei, the axial vector contribution is proportional to the nucleus spin [18], therefore it is negligible as compared to vector contribution.The cross section for scattering on the nucleus is: where F (q 2 ) is the nuclear form factor, , and m N is the nucleon mass.
The form factor, F (q 2 ), is a source of uncertainty in the scattering cross section, in particular at the highest nuclear recoil energy.Since there are no direct experimental mea-surements of the neutron distribution in the nucleus for CsI, F (q 2 ) may be estimated from the CEνNS data itself.Standard parameterizations are the Helm model [19], or the symmetrized Fermi model [20].Each of these model paramerizations are functions of the size of the neutron distribution, R n , which can be viewed as a parameter in each model.For each of these form factor parameterizations and including only SM contributions to the scattering rate, the best-fitting neutron radius from the COHERENT data is R n 5.5 fm [11].
The SM cross-section above may be modified due to a new mediating particle which couples to neutrinos and either electrons or quarks.For example, consider a new mediator Z µ with the following interaction terms: where g ν , g f,v , and g f,a are constants associated with new physics.With this interaction we can redefine the couplings (g v , g a ) of Eq. ( 1) in the following way where q 2 is the momentum transfer, M Z is the mass of the mediator, and the (−) sign applies for the case of anti-neutrino scattering.
As an extension, we can consider the case where the quarks couple to Z via a loop containing hidden sector particles [10].In this case a form factor parameterizes the quark coupling to the Z .This form factor scales as q 2 /Λ 2 , where Λ is the scale in the hidden sector associated with the mass of the mediator particle that generates the quark-hidden sector interactions.The form factor goes to unity when the quark couples directly to the Z .In this scenario, g ν g f in Eq. ( 4) is modified to g ν g f q 2 /Λ 2 , where Λ is the mass scale of the hidden sector.
Motivated by the flavor content of the COHERENT stopped pion beam, in this paper we focus on g e and g µ , associated, respectively, with electron and muon-type neutrino couplings with the Z .The up and down quark couplings to Z are represented by g u and g d , respectively, and we will assume that g u = g d = g ν = g.Under this assumption the NSI in the hidden sector scenario.Including this parameter will have the effect of reducing the predicted number of events as it interferes destructively with the SM weak charge.On the other hand, if we set g u = g d = −g ν = g, the NSI parameter will become large and negative.

III. DATA ANALYSIS
The COHERENT collaboration has published data on the energy and timing distribution of its nuclear recoil events [21].Figure 1 shows the energy and timing distribution of the events, with a classification of prompt and delayed events coming from the COHERENT collaboration.Note that the energy distribution is given in units of photoelectrons (PEs).
Our goal is to perform a joint analysis using both the timing and energy distribution to identify possible contributions from BSM physics.To represent BSM physics, we take as parameters the couplings g e , g µ , and mediator mass M Z .To represent the uncertainty in the nuclear structure, we take R n as a parameter and use the Helm form factor model. (right) as measured by COHERENT.The vertical axis in the time distribution gives the probability that an event is in a given time bin, while in the energy distribution it is the number of counts in that bin.In the energy distribution, the signal and background events are shown.
We define L as the likelihood function of the data given the model parameters.To form the likelihood function, we assume that the sum of the observed nuclear recoil plus background counts, N obs (t, E), at time t and energy E follows a Poisson model with parameter where In addition to the counts from the signal and the background components, Eq. ( 5) involves the uncertainty parameter α to account for the systematic uncertainties from flux, form factor, quenching factor, and signal acceptance uncertainties.Motivated by the results reported from COHERENT [1], we assume this parameter follows a normal distribution with zero mean and standard deviation σ α = 0.28.Defining θ = (g e , g µ , m, R n ) as the model parameters, the most general likelihood function for the parameters given the data is In comparison to previous analyses of the COHERENT data, this likelihood explicitly includes information from the timing distribution of the data.We take the bin width in energy space to be given by the bin width in the space of the number of photoelectrons, n e , then convert this to recoil energy space using the relation n e = 1.17(E/keV).For the timing data we take the bin width directly from the COHERENT data, ∼ 0.5µs.We obtain the posterior probability densities for the model parameters using the multinest package [22] with flat prior distributions on the parameters.
To test the robustness of our modeling in light of the steady-state (SS) background observed by COHERENT, we consider three separate scenarios for it.In the first scenario, which we refer to as background model (a), the SS background is fixed (N bg (t, E) = N bg,obs ) by the reported COHERENT measurements, and the model parameters are only (g e , g µ , M Z , R n ), and we do not need to integrate over the background in Eq. ( 6).In the second scenario, which we refer to as background model (b), we take the background shape from the COHERENT data, but allow for an overall scaling of the background to account for an uncertainty in the background normalization.Instead of integrating over the background for each bin, we only integrate over the total number of background counts in this scenario.In the third scenario, which we refer to as background model (c), we take the background to be a poisson model in each bin, and integrate over the predicted number of events in each bin according to Eq. ( 6).In this third scenario, in order to mitigate a bias in our reconstructed parameters, we exclude bins in which there are zero background counts.
To compare any two models, say models 0 and 1, defined by different sets of parameters, we consider the test statistic Here we define log(L 0 ) as the log-likelihood for model 0 in which only R n is free and other three parameters are set to zero, and log{L 1 ( θ)} is the log-likelihood for the model in which at least one of the parameters in (g e , g µ , M Z ) is free.For the latter model, θ denotes the maximum likelihood estimator (MLE) of θ.The p-value of the test is p = pr(χ 2 η > U obs ), where U obs is the observed value of U and χ 2 η is the chi-square distribution with η degree of freedom.Here, η is the difference between the number of estimated parameters in models 0 and 1.A small p-value provides significant evidence against the null hypothesis (SM).Applying this general procedure we can test if a model parameter (i.e., g µ or g e ) is positive.We define the corresponding significance Z as Φ −1 (1 − p/2), where Φ −1 is the inverse cumulative distribution function of the standard normal distribution.Significant results can also be seen through a large value of Z.

IV. RESULTS
We begin by considering how well R n is determined for each of our three different background models.The posterior probability distributions for R n are shown in the left panel of Figure 2 for each of these background models.The value of R n 5.5 fm obtained for our background model (c) is consistent with the result obtained in Ref. [11], though these authors obtained this result from background subtracted data.For background models (a) and (b), the best-fit R n is slightly lower than that from model (c), though is still statistically consistent with the result from model (c).
With the allowed parameter space for R n now understood, we move on to adding BSM physics in the form of light mediators.In order to simplify our analysis, for our BSM scenario we consider models with couplings such that g ν = g u = g d = g .To compare the sensitivity of the data to different flavors, we either fit for g µ and fix g e , or vice versa.The key features of our analysis are unchanged if we use a different relation among the couplings.The allowed parameter space shifts little bit towards larger values of R n if we consider the form factor ∼ q 2 /Λ 2 in the quarks-Z coupling.
In Figures 3 and 4, we show the resulting posterior probability distributions for g µ (as- suming that g e = 0) and similarly for g e (assuming that g µ = 0), for the cases of M Z = 10 and 1000 MeV.We show the result for background model (c), in which the background is taken to be a poisson model in each energy and time bin.The figures contain distributions using energy information alone, and distributions using both energy and timing information.For both M Z = 10 and 1000 MeV, the g µ distributions are better constrained when including timing data, and deviate from the SM prediction that these couplings are zero.
Though we have assumed background model (c) in Figures 3 and 4, we find that these distributions are relatively insensitive to the assumed background model.To quantify this, in Table I we show the significance Z for each of the three different background assumptions.
In Table I we also show the corresponding values of Z for an analysis with the energy data alone.Here we find that statistical significance of the deviation from the SM is lower for all background models, as compared to the analysis that utilizes both energy and timing data.The results in Table I show that the timing data provides additional information on the flavor content of the fluxes that is not provided by the energy data alone.The bin-0 .0 0 0 0 0 0 0 .0 0 0 0 2 5 0 .0 0 0 0 5 0 0 .0 0 0 0 7 5 0 .0 0 0 1 0 0 0 .0 0 0 1 2 5 0 .0 0 0 1 5 0 0 .0 0 0 1 7 5 0 .0 0 0 2 0 0 g by-bin likelihood analysis that we employ is able to statistically separate the prompt and delayed distributions, with the timing information more strongly constraining the prompt g µ component.In contrast, the g e component only contributes to delayed neutrino recoil spectrum, and is less well constrained when adding in the timing information.TABLE I: Significance in the likelihood ratio test for different assumed background models, as defined in the text.The first column gives the mass of the mediator that is assumed, with the first row representing the case in which the mediator mass is a free parameter.In all cases we take g e = 0 and g µ = 0.Each primary entry is obtained using the combined energy and timing data, and the numbers in parenthesis use just the energy data alone.
From the figure we can see for mediator mass from 1 MeV to 1 GeV the energy and timing data imply a best-fit deviation from the SM, while the energy data alone is more consistent with the SM.The shape of the boundary in both plots can be understood as follows: in the small mediator mass region, q 2 M 2 Z , the NSI parameter is independent of the small mediator mass, while in the large mediator mass region, M 2 Z q 2 , the NSI parameter depends on g 2 /M 2 Z , thus, in log space the slope is about 1.The isolated island at large mediator mass region is because the global degeneracy for the weak charge across all energy bins (since the NSI parameter is independent of energy).On the contrary, if a hidden sector is introduced to generate a form factor ∼ q 2 /Λ 2 , the NSI parameter becomes independent of energy in the smaller mass region and consequently, the degeneracy appears in the smaller mass region.
The above results favoring g = 0 are independent of the uncertainty in the assumed value for R n , since we have taken R n as a parameter in our analysis.It is however interesting to determine the manner in which R n is degenerate with the BSM model.As an example we can compare to the case in which g u = g d = −g ν = g, as in this case there is destructive interference with the SM weak charge.The resulting R N distribution with this assumption is shown in the right-hand panel of Fig. 2.Here we see that the degeneracy is such that larger R N is now more favored, though the results are statistically consistent with the assumption of g u = g d = g ν = g.We stress again that this distribution of R n is likely maximally conservative, and that input from the nuclear model of CsI will help break this degeneracy.

V. DISCUSSION AND CONCLUSION
In this paper we have performed a fit to the energy and timing distribution of nuclear recoil events from the COHERENT data.A combined timing and energy analysis provides information on the flavor content of the neutrino flux beyond what is obtained with the energy data alone.We have shown that including the information in both the energy and timing distributions of the COHERENT data, there is a ∼ 2σ deviation between the bestfitting model and the SM prediction.Light mediators in the mass range ∼ 10 − 1000 MeV are able to provide a good fit to the data.Though our primary analysis is performed within the context of models with g ν = g u = g d = g for light mediators, in the appropriate limits of large mediator mass, > ∼ 1 GeV, this analysis is equivalent to a heavy mediator NSI analysis that has been performed by the COHERENT collaboration and by previous authors.We also analyzed a BSM scenario where an additional factor ∼ q 2 /Λ 2 exists in the quarks-Z coupling arising from the hidden sector interactions.
The analysis that we have performed in this paper highlights that there is additional statistical information in an event-by-event analysis of COHERENT nuclear recoil events.
Given information on individual nuclear recoil events, our analysis may be extended in the future by considering an unbinned limit to the likelihood, and may include appropriate uncertainty on the energy and time of an event.Additional BSM physics may then ultimately be tested, e.g. in the form of sterile neutrinos or neutrino magnetic moment.
Future data from COHERENT, as well as from reactor-based experiments, will better clarify the results that we have obtained.An important expected near-term result is CO-HERENT data using different nuclear targets [23].In addition to increasing the exposure, different nuclear targets will help further decrease the uncertainties that arise due to the nuclear structure.

FIG. 1 :
FIG.1: Timing distribution for the prompt and delayed components (left) and energy distribution

2 FIG. 2 :
FIG. 2: Posterior probability distributions for the effective size of the nucleus R n .The left panel compares different background models, and the right panel compare different scenarios for new interactions.Here we assume the bin-by-bin fluctuating background model (c) as described in the text.The (blue) SM distribution only assumes R n as a parameter, while the BSM distribution takes (R n , g , and M Z ) as free parameters.

FIG. 3 :
FIG.3: Posterior probability distributions for g µ (top row) or g µ /Λ (bottom row, if there is form factor), using the energy data alone (orange) and using the combined energy and timing data (blue).The left column assumes a mediator mass of M Z = 10 MeV, and the right column assumes a mediator mass of M Z = 1000 MeV.

FIG. 5 :
FIG.5: Heat map of the probability density in log 10 (M Z ) vs log 10 g µ (top row) and log 10 (M Z ) vs log 10 (g µ /Λ) (bottom row, when there exists a form factor q 2 /Λ 2 ) parameter space using energy and timing data (left) and energy data alone (right).Here we assume background model (c).
Poisson model with parameter N bg (t, E).Moreover, in the absence of any prior information about N bg (t, E), we use a non-informative prior on it, so that π(N bg , E, ) is the number of neutrino-induced nuclear recoil events predicted from the theory and N bg (t, E) denotes the true background count.Note that by definition N bg (t, E) is not observed.Rather, we have observed background counts, denoted by N obs,bg (t, E), that are proxy for the true background counts.We assume that given N bg (t, E), N obs,bg (t, E) follows a