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

TL;DR

This paper combines energy and timing data from the COHERENT experiment to search for light mediators beyond the Standard Model, finding a mild preference for a 10-1000 MeV mediator with specific coupling strengths.

Contribution

It introduces a combined analysis of energy and timing data to constrain light BSM mediators, accounting for various uncertainties and highlighting the added value of timing information.

Findings

Data favor a 10-1000 MeV mediator over the Standard Model at < 2σ

Timing data provides unique flavor distribution information

Results are robust against uncertainties in nuclear and experimental parameters

Abstract

We search for beyond the Standard Model (BSM) physics by combining COHERENT energy and timing data. Focusing on light, $\lesssim$ GeV mediators, we find the data favor a $\sim 10-1000$ MeV mediator as compared to the Standard Model (SM) best fit at the $\lesssim 2\sigma$ level. The best-fit coupling range is $g \sim 10^{-5}-10^{-3}$. The timing data provides statistical information on the neutrino flavor distributions that is not attainable from the nuclear recoil energies alone. This result accounts for uncertainty in the effective size of the neutron distribution, and highlights the power of including uncertainties on experimental backgrounds, nuclear structure, and BSM physics.

Figures

Tables

Table 1. Table 1: 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 subscript 𝑔 𝑒 0 g_{e}=0 and g μ ≠ 0 subscript 𝑔 𝜇 0 g_{\mu}\neq 0 . Each primary entry is obtained using the combined energy and timing data, and the numbers in parenthesis use just the energy data alone.

Mediator mass, $M_{Z^{\prime }}$ (MeV)	Fixed (model (a))	Fixed shape (model (b))	Varying (model (c))
free	$1.4\left(0.7\right)$	$0.9\left(0.6\right)$	$1.1\left(0.6\right)$
$10$	$1.9\left(1.2\right)$	$1.4\left(1.1\right)$	$1.6\left(1.0\right)$
$100$	$1.9\left(1.1\right)$	$1.4\left(1.1\right)$	$1.6\left(1.1\right)$
$1000$	$1.9\left(1.2\right)$	$1.4\left(1.1\right)$	$1.6\left(1.1\right)$

Table 2. Table 2: The 1 ​ σ 1 𝜎 1\sigma constraining windows on g μ subscript 𝑔 𝜇 g_{\mu} and g e subscript 𝑔 𝑒 g_{e} using energy plus timing information. The results are shown using background model (c).

$M_{Z^{\prime }}$ (MeV)	10	100	1000
$g_{\mu }$	$[1.87,6.65]\times {10}^{-5}$	$[0.41,1.47]\times {10}^{-4}\oplus [2.47,2.66]\times {10}^{-4}$	$[0.48,1.32]\times {10}^{-3}\oplus [2.17,2.47]\times {10}^{-3}$
$g_e$	$[0,6.12]\times {10}^{-5}$	$[0,1.53]\times {10}^{-4}\oplus [2.53,2.84]\times {10}^{-4}$	$[0,1.22]\times {10}^{-3}\oplus [2.22,2.77]\times {10}^{-3}$

Equations

$$\frac{d\sigma }{dE}=\frac{G_F^{2}m}{2\pi }\left((g_v+g_a{)}^2+(g_v-g_a{)}^2{\left(1-\frac{E}{E_{\nu }}\right)}^2+(g_a^2-g_v^2)\frac{mE}{E_{\nu }^2}\right),$$

$$\text{derivative}\sigma E=\frac{G_F^2{Q}_V^2}{2\pi }{m}_N\left(1-\left(\frac{m_{N}E}{E_{\nu }^2}\right)+{\left(1-\frac{E}{E_{\nu }}\right)}^2\right)F(q^2)$$

$$\mathcal{L}\supset Z_{\mu }^{\prime }(g_{\nu }^{\prime }{\bar{\nu }}_L{\gamma }^{\mu }{\nu }_L+g_{f,v}^{\prime }\bar{f}{\gamma }^{\mu }f+g_{f,a}^{\prime }\bar{f}{\gamma }^{\mu }{\gamma }^{5}f),$$

$$(g_v,g_a)\Rightarrow (g_v,g_a)+\frac{g_{\nu }^{\prime }(g_{f,v}^{\prime },\pm g_{f,a}^{\prime })}{2\sqrt{2}{G}_F(q^2+M_{Z^{\prime }}^2)},$$

$$\lambda (t,E)=(1+\alpha )N(t,E,ϵ)+N_{bg}(t,E),$$

$$\mathcal{L}(\vec{\theta }|t,E)$$

Full text

I Introduction

The COHERENT collaboration has reported the first detection of coherent neutrino-nucleus elastic scattering (CE$\nu$NS) Akimov et al. (2017). 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, $\nu_{\mu}$, thus arrive from a prompt decay of charged pions, whereas $\bar{\nu}_{\mu}$ and $\nu_{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\pm 22$ nuclear recoil events from CE$\nu$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) Ohlsson (2013); Miranda and Nunokawa (2015) which may arise from heavy or light mediators Coloma et al. (2017a, b); Liao and Marfatia (2017); Dent et al. (2018); Billard et al. (2018); Lindner et al. (2017); Farzan et al. (2018); Brdar et al. (2018). The data also provide novel constraints on new physics that manifests through the neutrino sector, including generalized scalar and vector neutrino interactions Aristizabal Sierra et al. (2018), as well as hidden sector models Datta et al. (2019). It also sets independent constraints on the effective neutron size distribution of CsI Ciuffoli et al. (2018); Aristizabal Sierra et al. (2019); Papoulias et al. (2019), and on sterile neutrinos Kosmas et al. (2017); Blanco et al. (2019). Low mass mediator models are particularly interesting since there are no Large Hadron Collider (LHC) constraints for mediator masses $\lesssim$ GeV, and these models connect to new ideas associated with sub-GeV dark matter in the range MeV to GeV Essig et al. (2013). 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 $\nu_{\mu}$ 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 Denton et al. (2018). 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 consider a scenario in which the quark-$Z^{\prime}$ 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 $\lesssim 2\sigma$ level. We also consider constraining the LMA-Dark solution using the timing data.

II Light mediators and Non-Standard Interactions

In this section we briefly discuss the SM prediction for CE$\nu$NS. We introduce a model for light mediators as an example of BSM physics, and highlight how this model alters the predicted CE$\nu$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

[TABLE]

where $G_{F}$ is the Fermi constant, $E$ is the recoil energy, $E_{\nu}$ is the incident neutrino energy, $(g_{v},g_{a})=(T_{3}-2Q_{\rm em}{\rm{sin}}^{2}\theta_{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_{\rm em}$ is the electromagnetic charge, and $\theta_{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\lesssim\mathcal{O}(50)$ MeV), coherent scattering occurs. For coherent elastic scattering off of nuclei, the axial vector contribution is proportional to the nucleus spin Barranco et al. (2005), therefore it is negligible as compared to vector contribution. The cross section for scattering on the nucleus is:

[TABLE]

where $F\left(q^{2}\right)$ is the nuclear form factor, $Q_{V}\equiv Z\left(2g^{u}_{v}+g^{d}_{v}\right)+N\left(g^{u}_{v}+2g^{d}_{v}\right)$, 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 measurements of the neutron distribution in the nucleus for CsI, $F(q^{2})$ may be estimated from the CE$\nu$NS data itself. Standard parameterizations are the Helm model Lewin and Smith (1996), or the symmetrized Fermi model Piekarewicz et al. (2016). 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}\simeq 5.5$ fm Ciuffoli et al. (2018).

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_{\mu}^{\prime}$ with the following interaction terms:

[TABLE]

where $g_{\nu}^{\prime}$, $g_{f,v}^{\prime}$, and $g_{f,a}^{\prime}$ 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

[TABLE]

where $q^{2}$ is the momentum transfer, $M_{Z}^{\prime}$ 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^{\prime}$ via a loop containing hidden sector particles Datta et al. (2019). In this case a form factor parameterizes the quark coupling to the $Z^{\prime}$. This form factor scales as $q^{2}/\Lambda^{2}$, where $\Lambda$ 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^{\prime}$. In this scenario, $g_{\nu}^{\prime}g_{f}^{\prime}$ in Eq. (4) is modified to $g_{\nu}^{\prime}g_{f}^{\prime}q^{2}/\Lambda^{2}$, where $\Lambda$ 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_{\mu}$, associated, respectively, with electron and muon-type neutrino couplings with the $Z^{\prime}$. The up and down quark couplings to $Z^{\prime}$ are represented by $g_{u}$ and $g_{d}$, respectively, and we will assume that $g_{u}=g_{d}=g_{\nu}=g$. Under this assumption the NSI parameter defined as $\epsilon=g_{\nu}^{\prime}g_{f,v}^{\prime}/\{2\sqrt{2}G_{F}\,(q^{2}+M_{Z^{\prime}}^{2})\}=g^{2}/\{2\sqrt{2}G_{F}\,(q^{2}+M_{Z^{\prime}}^{2})\}$ is positive, or $\epsilon=g^{2}q^{2}/\{2\sqrt{2}G_{F}\,(q^{2}+M_{Z^{\prime}}^{2})\Lambda^{2}\}$ 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_{\nu}=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 Akimov et al. (2018a). 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_{\mu}$, and mediator mass $M_{Z^{\prime}}$. To represent the uncertainty in the nuclear structure, we take $R_{n}$ as a parameter and use the Helm form factor model.

We define $\mathcal{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

[TABLE]

where $N(t,E,\epsilon)$ 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}\left(t,E\right)$ 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}\left(t,E\right)$, $N_{obs,bg}(t,E)$ follows a Poisson model with parameter $N_{bg}\left(t,E\right)$. Moreover, in the absence of any prior information about $N_{bg}(t,E)$, we use a non-informative prior on it, so that $\pi(N_{bg}(t,E))\propto 1$ for $N_{bg}(t,E)\in[0,\infty)$.

In addition to the counts from the signal and the background components, Eq. (5) involves the uncertainty parameter $\alpha$ to account for the systematic uncertainties from flux, form factor, quenching factor, and signal acceptance uncertainties. Motivated by the results reported from COHERENT Akimov et al. (2017), we assume this parameter follows a normal distribution with zero mean and standard deviation $\sigma_{\alpha}=0.28$. Defining $\vec{\theta}=(g_{e},g_{\mu},m,R_{n})$ as the model parameters, the most general likelihood function for the parameters given the data is

[TABLE]

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/\textrm{keV})$. For the timing data we take the bin width directly from the COHERENT data, $\sim 0.5\mu$s. We obtain the posterior probability densities for the model parameters using the multinest package Feroz et al. (2009) 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_{\mu},M_{Z^{\prime}},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 $U=-2[\log(\mathcal{L}_{0})-\log\{\mathcal{L}_{1}(\hat{\theta})\}]$. Here we define $\log(\mathcal{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\{\mathcal{L}_{1}(\hat{\theta})\}$ is the log-likelihood for the model in which at least one of the parameters in $(g_{e},g_{\mu},M_{Z^{\prime}})$ is free. For the latter model, $\hat{\theta}$ denotes the maximum likelihood estimator (MLE) of $\vec{\theta}$. The $p$-value of the test is $p=\hbox{pr}(\chi^{2}_{\eta}>U_{\rm obs})$, where $U_{\rm obs}$ is the observed value of $U$ and $\chi^{2}_{\eta}$ is the chi-square distribution with $\eta$ degree of freedom. Here, $\eta$ 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_{\mu}$ or $g_{e}$) is positive. We define the corresponding significance $Z$ as $\Phi^{-1}(1-p/2)$, where $\Phi^{-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}\simeq 5.5$ fm obtained for our background model (c) is consistent with the result obtained in Ref. Ciuffoli et al. (2018), 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_{\nu}=g_{u}=g_{d}=g^{\prime}$. To compare the sensitivity of the data to different flavors, we either fit for $g_{\mu}$ 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 $\sim q^{2}/\Lambda^{2}$ in the quarks-$Z^{\prime}$ coupling.

In Figures 3 and 4, we show the resulting posterior probability distributions for $g_{\mu}$ (assuming that $g_{e}=0$) and similarly for $g_{e}$ (assuming that $g_{\mu}=0$), for the cases of $M_{Z^{\prime}}=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^{\prime}}=10$ and $1000$ MeV, the $g_{\mu}$ 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 1 we show the significance $Z$ for each of the three different background assumptions. In Table 1 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 1 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-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_{\mu}$ 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.

Figure 5 shows the probability density in $\log_{10}(M_{Z^{\prime}})$ vs $\log_{10}(g_{\mu}$) or $\log_{10}(g_{\mu}/\Lambda)$ space. 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}\gg M_{Z^{\prime}}^{2}$, the NSI parameter is independent of the small mediator mass, while in the large mediator mass region, $M_{Z^{\prime}}^{2}\gg q^{2}$, the NSI parameter depends on $g^{2}/M_{Z^{\prime}}^{2}$, 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 $\sim q^{2}/\Lambda^{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\neq 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_{\nu}=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_{\nu}=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.

To quantify the constraints on the NSI parameters in Fig. 3 and Fig. 4, we show in Table 2 the $1\sigma$ window of $g_{\mu}$ and $g_{e}$ assuming 10, 100, and 1000 MeV mediator masses using our results from the energy and timing analysis. Note that this range is consistent with general constraints on the couplings and masses of light mediators Bauer et al. (2018).

Constraints on the NSI parameters also come from neutrino oscillation experiments. To obtain constraints from neutrino oscillation experiments, the NSI parameters are expressed as $\epsilon=g^{2}/\{2\sqrt{2}G_{F}\,(M_{Z^{\prime}}^{2})\}$. Since the neutrino undergoes forward scattering off matter, it is independent of $q^{2}$. The existing constraints from the neutrino oscillation experiments are consistent with our result, within the $\sim 2\sigma$ range Coloma et al. (2017a). For example, the allowed range for a heavy mediator for $\epsilon^{u}_{\mu\mu}$ is $[0,0.005]$ (our model only has positive $\epsilon$). This converts to $g=[0,0.0013]$ at $M_{Z^{\prime}}=1~{}\mathrm{GeV}$, which covers the $1\sigma$ region as shown in Fig. 3.

There also exists so-called a LMA-Dark solution, which is characterized by larger NSI couplings and is not strongly constrained by oscillation experiments. Ref. Denton et al. (2018) discusses this scenario specifically in the context of the COHERENT data. Since the SM is not nested within the LMA-Dark solution, to establish significance in this case we must appeal to the Bayes factor. We find that the Bayes factor is $1.12~{}(0.85)$ with (without) timing information, which means this kind of NSI is also slightly preferred over SM if we include the timing data.

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 $\sim 2\sigma$ deviation between the best-fitting model and the SM prediction. Light mediators in the mass range $\sim 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_{\nu}=g_{u}=g_{d}=g^{\prime}$ for light mediators, in the appropriate limits of large mediator mass, $\gtrsim 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 analyzed a BSM scenario where an additional factor $\sim q^{2}/\Lambda^{2}$ exists in the quarks-$Z^{\prime}$ coupling arising from the hidden sector interactions. We showed the $1\sigma$ allowed window of $g$ in our model. Our result is consistent with dark photon searching experiments as well as oscillation experiments. We also showed the energy plus timing information can be used to constrain LMA-Dark solution of neutrino oscillation, and that LMA-Dark solution is also prefered than SM with Bayes factor $1.12$.

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 COHERENT data using different nuclear targets Akimov et al. (2018b). In addition to increasing the exposure, different nuclear targets will help further decrease the uncertainties that arise due to the nuclear structure.

ACKNOWLEDGMENTS

BD and LES acknowledge support from DOE Grant de-sc0010813. We acknowledge support from COS-STRP (TAMU), and would like to thank Grayson Rich and Kate Scholberg for discussions on this paper.