Probing $L_\mu-L_\tau$ models with CE$\nu$NS: A new look at the combined COHERENT CsI and Ar data

The minimal gauged $U(1)_{L_\mu-L_\tau}$ model has long been known to be able to explain the tension between the theoretical and experimental values of the muon magnetic moment. It has been explored and tested extensively, pushing the viable parameter space into a very tight corner. Further, embedding the $U(1)_{L_\mu-L_\tau}$ model in a supersymmetric (SUSY) framework has been shown to relax some of these constraints and has recently been shown to explain the electron anomalous magnetic moment as well. In this model, the logarithm of the mass ratio of third to second generation (s)leptons control the non-negligible kinetic mixing and may crucially alter many of the constraints. We confront both the non-SUSY and SUSY versions of this class of models with the CsI(2017), the recently released CENNS10 data from the liquid Argon detector as well as the updated CsI(2020) data of the COHERENT experiment. We use the recoil energy and timing binned data from CsI(2017) and the energy, time, and Pulse Shape Discriminator binned data from CENNS10 to find estimates for the model parameters in a likelihood maximization test. We also show updated exclusions using all of the above data from the COHERENT Collaboration, as well as projected exclusions from the ongoing Coherent CAPTAIN-Mills experiment. The $(g-2)_\mu$ favored values of the $U(1)_{L_\mu-L_\tau}$ gauge coupling that are still unconstrained overlap with the estimates from COHERENT data within $1\sigma$. The combined COHERENT data is found to prefer the presence of the $U(1)_{L_\mu-L_\tau}$ gauge boson over the Standard Model at $\sim1.4\sigma$. The global minima of a chi-square deviation function using CsI(2020) as well as CENNS10 total counts has significant overlap with the $(g-2)_{\mu}$ favored parameter space in the context of the SUSY and non-SUSY $L_{\mu}-L_{\tau}$ models with a mediator mass in the $20-100$ MeV range.


I. INTRODUCTION
The prospect of observing coherent elastic neutrinonucleus scattering (CEνNS) was first proposed roughly 47 years ago [1]. However, the crucial discovery was of the fact that the scattering cross-section involved with the CEνNS process would be many times larger than the corresponding process with a single nucleon. The CEνNS cross-section was expected to scale approximately with the square of the neutron number of the nuclei [1,2].
In spite of the large predicted event rates, which meant that detectors with masses as small as a few kilograms could be used to detect it, CEνNS evaded observation for four decades, primarily because of the small (few keV) nuclear recoil energies involved. The COHER-ENT Collaboration finally reported its observation of the process with stopped-pion neutrinos on a CsI detector, 43 years after it was first proposed, and the recoil energy distribution agreed with the Standard Model (SM) within 1σ. The Collaboration also reported that the twodimensional (energy and time) profile, fit to the maximum likelihood, prefers the presence of CEνNS against its absence at the 6.7σ confidence level [3]. Henceforth, we refer to this as the CsI(2017) dataset. This was followed up by the observation of CEνNS on a liquid Ar detector which also rejects the absence of CEνNS at 3.5σ with 3D-binned data in their recoil energy, timing, as well as pulse-shape-discriminator (PSD) data [4]. The CO-HERENT Collaboration has also recently announced the impending publication of an update to the initial dataset from their CsI detector [5], referred to as CsI(2020), which rejects the absence of CEνNS at 11.6σ. Not only have these results added to the understanding of neutrino interactions and nuclear form factors [6][7][8][9][10][11][12] within the SM, but they have also provided a benchmark to test any propositions of physics beyond the SM that altered them.
The modification to the SM CEνNS rate due to the L µ −L τ gauge boson proceeds via the kinetic mixing ( ), a necessarily small parameter ( 10 −5 ). As a consequence, the modification is dominated by the interference term, which is destructive as long as the sum of the u-and d-quark beyond the Standard Model (BSM) couplings to the neutrinos remain positive. This holds for all but a few pathological parameter choices for the U (1) Lµ−Lτ model. This is owing to the fact that these couplings in the L µ − L τ model are proportional to Q em . While the presence of the L µ − L τ gauge boson thus serves to reduce the CEνNS rate, the CsI detector has observed a count lower than that predicted by the SM and the Ar detector has reported an opposite trend. It is true that both results are within 1σ of the SM prediction, still, a more rigorous analysis involving the full energy as well as timing data has shown that the CsI data prefers new physics scenarios that lead to a reduction of the number of counts at ∼ 2σ [34]. We re-evaluate the situation in the light of the recent Ar data adding a pull in the opposite direction.
Our analysis finds that although the preference for destructive BSM effects has a marked reduction, it still remains slightly larger than 1σ. We use our analysis to show updated exclusions on the L µ −L τ parameter space. The exclusions from the Ar data are much stronger than the corresponding exclusions from CsI(2017), however the exclusions from the combined data adhere closely to those from the CsI(2017) data. The CsI(2020) combined with the Ar data lays much stronger bounds than the former, still; they are not the strongest constraints in the relevant parameter space. However, the preferred region of parameter space from the likelihood test, as well as the global minima of the chi-square deviation has significant overlap with the (g − 2) µ preferred region in the context of the L µ − L τ model. In addition, we also show projected exclusions from the Coherent CAPTAIN-Mills (CCM) [81,82] experiment which is another experiment looking for CEνNS with stopped-pion neutrinos on a liquid Ar detector. In comparison with the COHERENT experiment, CCM uses a less intense neutrino beam but has a much larger detector fiducial volume (∼5 tons as opposed to ∼24 kgs at COHERENT). In case the CCM observations correspond to the SM predictions, the exclusions would be the strongest to date in the pertinent region of the parameter space.

II. CEνNS AT COHERENT
The COHERENT Collaboration was designed to detect and study CEνNS over a range of detectors. It uses a high intensity neutrino beam produced at the Spalla- FIG. 2. CEνNS energy (2(a)), timing (2(b)) and F90 (2(c)) distributions from the CENNS10 detector [4]. Also shown in Fig.2(a) is our SM CEνNS prediction for CENNS10 along with the associated total uncertainty shown shaded green.
tion Neutron Source (SNS) of the Oak Ridge National Laboratory. A 1 GeV proton beam incident upon a mercury target at 60 Hz produces ∼ 0.08 − 0.09 neutrinos per flavor per proton on target (POT). Upon hitting the Hg target, the proton beam produces pions. The negative pions are captured almost entirely within the target while the positive pions decay to produce the neutrinos, The monoenergetic ν µ from the first pion decay are "Prompt" and arrive at the target within a very short span of time ( 1.5µs) after the POT. The ν e andν µ from the subsequent muon decay are "Delayed" and have an energy profile as they are produced in a three-body decay. The neutrino flux from the SNS may be described by, where r = the number of neutrinos per flavor produced per proton on target, N POT = the total number of protons on target during the data-taking period, and L = the distance of the detector from the neutrino source. This results in a distinct timing spectra for the incoming neutrinos at the SNS as shown in Fig.1(b). Data from the observation of CEνNS at two different detectors has been made public by the COHERENT Collaboration to date: the CsI and the CENNS10 liquid Ar detectors.

A. Backgrounds
For both detectors, there are two broad classes of backgrounds that affect the analysis. The first, steady state backgrounds, are present in both beam-on and beamoff data and may be handled using prescribed methods of background subtraction as described in Refs. [4,35]. The second, beam-related-neutrons (BRN) and neutrinoinduced-neutrons (NIN) occur in time with the SNS beam and are hence more difficult to separate. The BRN can be further classified into prompt and delayed BRN, the former are the fast neutrons produced by the SNS target itself while the latter are produced by the neutrinos interacting with the detector shielding. Both of these backgrounds generate events similar in nature to the CEνNS events and need to be correctly accounted for. The COHERENT Collaboration uses the different recoil energy and timing distributions for these events to discriminate between this background and the CEνNS signal.
The COHERENT Collaboration has released data from dedicated measurements of both the kinds of backgrounds for CsI as well as CENNS10 detectors. In our analyses we compare the total observed counts with the model prediction added to these background counts using appropriate uncertainties for all the distributions. This is consistent with the method used by the Collaboration itself to generate the best-fit CEνNS and background counts. As a cross-check, we use our method to re-evaluate these best-fit counts and find a very close match with those reported by the Collaboration. efficiency of detection and background distinction was poor below ∼ 10 keV. It was in place from 2015 to 2019 and was used to make the first observation of CEνNS with 171.7 beam-on live-days, which amounted to approximately 1.76×10 23 protons on target. The quenching factor for the CsI detector, describes the number of photoelectrons detected by the photo multiplier tubes per keV nuclear recoil energy produced in a CEνNS event. Reference [35] describes the detector acceptance efficiency in terms of the detected number of photoelectrons (x) by the function where, a = 0.6655 +0.0212 Once the number of photoelectrons is mapped to nuclear recoil energies using Eq.3, we get the acceptance A(E r ). The CsI detector reported [3] an observation of 134 ± 22 CEνNS events at a 6.7σ C.L. against an SM prediction of 173 ± 48 events in 2017. The energy, as well as the timing distribution [the probability distribution function (PDF) generated from observation] of the observed data, as reported, is shown in Fig.1. This observation has a total fractional uncertainty of 28% which is a combination of neutrino flux (10%), nuclear form factor (5%), quenching factor (25%), and cut-acceptance (5%) uncertainties. The COHERENT Collaboration has recently announced the publication of additional data from the CsI detector [5], however only the total CEνNS and background counts have been made public along with the respective uncertainties. We use the public data from the first COHERENT data release in the likelihood based analysis (discussed in Sec.V) and the recently announced total counts in a χ 2 test.

C. CENNS10
This is a liquid Argon detector with 24Kg of active fiducial mass placed 27.5m away from the SNS target off the proton beam axis by ∼ 135 • and operational since 2017. With a threshold nuclear recoil energy of ∼ 3 keV, the detector recorded data worth approximately 1.38 × 10 23 POT. Reference [4] describes the quenching factor, The acceptance, A(E r ), for the CENNS10 detector is interpolated from efficiency data binned in nuclear recoil energy provided with the data release. In addition to the energy and timing information of the observation, this data release also includes the PSD distribution. This observable is used to distinguish between electron and nuclear recoils in Ar. More on the PSD variable (F90) and the cuts used on it can be found in Ref. [4]. We use F90 cuts identical to Analysis A described in this reference.
The CENNS10 detector reported [4] the observation of 159 ± 43(stat) ± 14(syst) CEνNS events at a significance of 3.5σ, as opposed to a SM prediction of 128 ± 17 events. The energy, timing and F90 distributions from the COHERENT Collaboration are shown in Fig.2. The involved uncertainties were divided into two part. The first set combines the quenching factor (1%), calibration (0.8%), detector efficiency (3.6%), prompt light fraction (7.8%), form factor (2%), and neutrino flux (10%) uncertainties for a total of 13%. This is an uncertainty on the CEνNS rate prediction and is called the prediction uncertainty for the rest of this work. The second set includes uncertainties in the F90 energy dependence(4.5%), neutrino arrival time(2.7%), BRN energy shape(5.8%), BRN arrival time-mean(1.3%), and BRN arrival timewidth(3.1%) and leads to a total uncertainty of 8.5% on the CEνNS best-fit value. This set will be referred to as the fit uncertainty. Public data from the CENNS10 detector is used in a likelihood as well as an χ 2 test in conjunction with the observations from the CsI detector.

D. Update to CsI data
The COHERENT Collaboration has recently announced the publication of additional data from its CsI detector [5] that has observed 306 ± 20 events against an SM prediction of 333±11(th)±42(ex). While the full details of the observation have not been made public yet, we have enough information to conduct a preliminary analysis with it. The Collaboration has updated its signal prediction as well as the CEνNS data analysis which has lead to a significant reduction in the overall uncertainty of the signal from 28% to 13%. This has also resulted in a reduction in both the predicted and observed CEνNS count.
The update releases data taken on the CsI[Na] detector until June 2019, which amounts to ∼ 2.2× more protons on target compared to the previous CsI data release. In the absence of detailed quenching factor and efficiency information from the new analysis, we use those from the previous data release. According to the Collaboration, the new analysis leads to a ∼ 9% reduction in the SM prediction for the CsI(2017) dataset. We find that the predicted signal for the CsI(2020) dataset using the earlier methods requires ∼ 14% reduction to match the new SM prediction released by the Collaboration. We also have information on the steady state background (1273 ± 38), BRN (17 ± 4), and NIN (5 ± 2), and the corresponding uncertainties for the new data.
Since we have only the total observed counts for both the signal and background in the new dataset, adopting a full likelihood maximization routine like the ones for the CENNS10 and old CsI data is overkill. Instead, we use a chi-square fit of the model parameters to the observed data in order to get the updated exclusions on the M Z − g X parameter space, with N obs = 306 denoting the total observed CEνNS count and N th denoting the model prediction. The quantities σ 2 obs = (δN obs ) 2 +(δN SS ) 2 +(δN BRN ) 2 +(δN N IN ) 2 and σ α = 0.13 parametrize the total uncertainties on the observed and predicted counts, respectively, δN obs = 20 In order to combine the CENNS10 data with this, we must define a similar chi-square function for it, with N obs = 159 and N th denoting the model prediction for the CENNS10 analysis. The quantities σ 2 obs = (δN obs ) 2 + (δN SS ) 2 + (δN pBRN ) 2 + (δN dBRN ) 2 and σ β = 0.1328 are the total uncertainties on the observed and predicted counts as for the previous case, In our analysis we define the total chi-square deviation as the sum of the two individual ones and minimize the chisquare deviation over the α and β nuisance parameters.

III. CEνNS AT CCM
The proposed CCM experiment is designed to detect CEνNS on a liquid Argon detector with a ∼ 5 ton fiducial mass using neutrinos from the Lujan target at the Los Alamos National Laboratory [81][82][83]. The detector is placed 20 m from the Lujan target that produces ∼ 0.04 neutrinos per flavor per POT. The Lujan facility uses an 800 MeV proton beam impinging upon a tungsten target at 20 Hz to produce pions that decay and produce neutrinos with the same energy shape as those at the SNS.
We assume 5000 hours of operation per year, for three years, for our analysis which is approximately 3 × 10 22 POT for the entire period. In the absence of any quenching factor, efficiency, or background information we only use a detector threshold of 25 keV and the CEνNS count alone for evaluating the projected exclusions. The detector efficiency above the threshold was kept at 100% and a 1 PE/keVnr quenching factor was used for our analysis. This was done for the sake of easy scalability once the full data from the detector becomes available.

IV. SIGNAL PREDICTION
In this section we shall describe the CEνNS cross section as it arises in the SM, discuss how it may be modified in the L µ − L τ model (in both the presence and absence of SUSY), and describe how we generate the signal prediction for the experiments.
In coherent neutrino-nucleus scattering events, the momentum transfer between the neutrino and the nucleus is smaller than or comparable to inverse the size of the nu- The axial vector contribution to this process is proportional to the nucleus spin [13], and hence is negligible in comparison to the vector contribution. The cross section for CEνNS scattering on a nucleus of mass m N , atomic number Z, and neutron number N is We use g v = T 3 − 2Q em sin 2 θ W , the quark vector coupling to the Z-boson where T 3 is the third component of the weak isospin, Q em is the electromagnetic charge, and θ W is the weak mixing angle.
The form factor F (q) is a description of the nucleon distribution within the nucleus and follows from its Fourier transform. It is captured by several parametrizations, e.g., the Helm [84], Symmetrized Fermi [85], and Klein Nystrand [86] approaches. Sticking to the conventions adopted by the COHERENT Collaboration in their own analysis of CsI and CENNS10 data, we use the Helm and Klein-Nystrand parametrizations for the two detectors, respectively. Note, however, that the choice of form factor parametrization has little or no effect on the results obtained.
The Helm form factor can be described by where q= √ 2m N E r is the absolute value of exchanged three-momentum, s=0.9 fm is the nuclear skin width and j 1 (x) = sin x/x 2 −cos x/x is the spherical Bessel function of order one. The neutron charge radius is given by Similarly, the Klein-Nystrand form factor is described by, The parameter a K = 0.7 fm enters the expression for the neutron charge radius, The size of the neutron distribution radius R n is taken to be 4.7 fm and 4.1 fm for the analyses involving CsI and Ar, respectively [8]. The SM CEνNS cross section may be modified in the presence of an additional gauge boson coupling to the neutrinos and/or quarks. The L µ − L τ gauge boson has direct couplings to second-and third-generation leptons as well as their scalar superpartners in the presence of SUSY. It also couples to any other charged field e.g. the quarks, via its kinetic mixing with the photon. It has the following interactions with the neutrinos and the quarks: Here Q X are the charges of the neutrinos under the additional gauged symmetry as defined in Ref. [80], and is the kinetic mixing parameter. It modifies the CEνNS cross section via the kinetic mixing of the Z with the photon [see Fig.3(a)]. In order to implement the modification of the SM CEνNS rate by the additional contributions it is enough to make the replacement in Eq. 10. This replacement modifies the quantity Q V in Eq. 10 which may be captured by the transformation and A, Z, and N are the mass number, atomic number, and neutron number, respectively, for the target nucleus (Cs, I, or Ar in our case, depending on the detector). The two terms in Eq.17 represent the u-quark and dquark couplings to the neutrinos multiplied by (A+Z) and (A+N) respectively and differ only by their electromagnetic charge. Note that in our case, owing to the inherent smallness of and hence Q BSM V , the interference term 2Q SM V Q BSM V in the differential cross-section dominates over the pure BSM part of it. Moreover, since Q SM V ∼ −0.5N , the interference is always destructive as long as the right hand side of eqn.17 remains positive.
This holds for the L µ − L τ model except for a few pathological parameter choices since the u-ν BSM coupling is positive and twice the d-ν coupling, which is negative. Thus, the CEνNS count suffers a reduction in this model. Figure 3(b) illustrates this effect for a choice of the L µ − L τ parameters. The black solid line in this figure denotes the SM prediction of the differential CEνNS event rate at the CENNS10 detector of the COHERENT experiment. The red dashed line shows the same in the presence of a 20 MeV L µ − L τ Z with the gauge coupling (g x ) set to 5.6 × 10 −4 and the SUSY contribution to taken to be vanishing. This choice of parameters can explain (g − 2) µ,e in the presence of SUSY [80] and does not predict an observable charged lepton flavor-violating decay width, even when neutrino masses and mixing may be explained in conjunction [78].
The presence of kinetic mixing that mediates the Z contribution to CEνNS may be assumed to be vanishing at the tree level, but it is unavoidably generated radiatively. In particular, the one-loop kinetic mixing in this class of models may be non-negligible in spite of superheavy BSM fields entering the loops [87] This expression for involves both charged leptons and their scalar superpartners. As discussed in Ref. [87], this quantity stands out in its dependence on only the ratio of the (s)lepton masses. While the ratio of the second and third generation leptons are fixed, we take the ratio r = mτ /mμ as an additional model parameter. In the absence of SUSY, only the first term contributes, as would be the case for degenerate sleptons (r = 1) even in its presence. It also assumes left-right degeneracy amongst the charged sleptons for the sake of simplicity. Relaxing these assumptions would only lead to the dependence of on additional mass ratios beyond r. We use Eq.18 to provide the value of in the replacement defined in Eq.16 in order to generate the modified differential cross section. For both Eqs.16 and 18, the transferred momentum is q 2 = 2M N E r , where M N is the nuclear mass and E r is the nuclear recoil energy. For the mediator mass (M Z ) regime that is pertinent for the L µ − L τ model ( 1 GeV), q 2 is comparable to M 2 Z . Unlike general neutrino NSIs, the BSM modification to the CEνNS rate is thus nuclear recoil energy dependent. The signal shape in the E r space due to this variation is taken into account by our likelihood function in the ensuing comparison. For a study of the impact of CEνNS observations on neutrino NSI parameter space see, for example, Ref. [88].
From the differential cross section, the total event rate in each recoil energy bin may be evaluated using, (19) where N D represents the number of target nuclei in the detector mass and is given by where, m det is the detector mass, N A is the Avogadro's number, m mol is the molar mass of the detector molecule and g mol is the number of atoms in a single detector molecule.
The presence of beyond-SM coupling of the neutrinos to quarks modifies the SM CEνNS rate. For opposite sign couplings of the neutrinos to up-and down-type quarks, as is the case in the presence of an L µ − L τ gauge boson, this effect is destructive as long as the two couplings are similar in magnitude. As a result, the presence of the L µ − L τ gauge boson serves to reduce the SM CEνNS rate for much of the parameter space, especially for low nuclear recoil energies [see Fig.3(b)]. Once the nuclear recoil energy spectrum has been generated, we use the timing and F90 distributions shown in Figs.1(b), 2(b), and 2(c) to generate the full experimental prediction for the respective detectors.

V. LIKELIHOOD TEST
This section describes the statistical analysis comparing the model predictions with experimental data from the COHERENT Collaboration. We adopt the method described in Ref. [34] for the analysis of the CsI data using both energy and time binning. We build upon this to devise a method to test our model against the full 3Dbinned CENNS10 data as well, the goal being an analysis using combined CsI and CENNS10 data binned in recoil energy and timing as well as PSD observables.
While defining the likelihood function for the CsI as well as CENNS10 data, we assume that the total count (CEνNS + background) follows a Poisson distribution in each bin. The model prediction added to the observed background counts acts as the parameter for the Poisson distributions used. The likelihood function is marginalized over the uncertainties in each of the separate distributions assuming a Gaussian envelope. In case of the CsI analysis, the observed background count is taken to follow a Poisson distribution with a hypothetical parameter N BG , that is marginalized. In the absence of clear uncertainties on the observed background count, this approach takes into account the statistical errors associated with this observation. Note that the uncertainties in each bin are marginalized over separately while building the likelihood function. This ensures that the signal shape information over all the free parameters ( θ) is accounted for while comparing the predicted and observed data.
The likelihood function for the CsI data given model where, are the parameter for the Poisson distribution, the Poisson distribution probability mass function (PMF) and the Gaussian distribution PMF respectively. N ( θ) is the CEνNS model prediction in each bin. The dummy variable α is used to marginalize the likelihood in each bin over the prediction uncertainty using a Gaussian envelope with standard deviation σ α = 0.28. The likelihood function for the CENNS10 data is defined with the same basic principles in mind. Unlike the CsI data, it has significantly large BRN observations that must be taken into consideration. The CENNS10 Datarelease [4] also provides the uncertainties associated with each of the three background components, that can now be incorporated directly. The likelihood function for the CENNS10 data is defined to be, Here, the parameter λ(α i , θ) is (24) where N SS , N pBRN , and N dBRN are the steady state, prompt BRN and delayed BRN background counts, respectively, in each bin. Just as in the case for CsI, the parameters, α i , are each the dummy variables that are marginalized to replicate uncertainties in each of the respective distributions as described in the first column of Table I. Only the prediction uncertainty (13%) is considered for the CEνNS count when estimating the model parameters.
In order to check the robustness of L Ar , we use it in repeating the likelihood maximization test described in Analysis A of Ref. [4]. Instead of marginalizing over the parameters α i , we pass them as parameters to be estimated by maximizing the likelihood function in Eq. 23 without the Gaussian envelopes. We use the SM prediction of CEνNS for the values of N ( θ) and a uniform prior on its normalization, α 1 . For the backgrounds, we use predicted counts instead of the best-fit counts as central values (as provided in Table I of Ref. [4]) and Gaussian priors on their normalizations with standard deviations The standard deviation for the prior on the steady state background normalization was taken to be the same as the uncertainty on its predicted value. The results of this table are shown in the second column of Table I. For comparison, we also show the corresponding bestfit values and the errors as reported by the COHERENT Collaboration and we find them to be an excellent match. Coming back to the estimation of model parameters using the likelihood functions defined in Eqs. 21 and 23, the combined likelihood function was defined to be the product of the two separate ones, We use this likelihood function in two sets of maximization tests, one where θ ≡ (g x ) with M Z as a fixed parameter and another where θ ≡ (g x , M Z ). The ratio r is taken as a fixed parameter for both analyses. We always use noninformative priors on the estimated parameters ( θ) and maximize the likelihood using the MultiNest algorithm [89]. We use the log-ratio test [90] to calculate the deviation of the maximum likelihood estimate (MLE) from the SM or the experimental best-fit observation while evaluating exclusions. In comparisons of two different models with different sets of estimated parameters, we define the test statistic, For our purposes, L 0 represents the SM likelihood with no estimated parameters and L(θ) is the likelihood at the MLE. According to Wilks' theorem this test statistic is distributed as a χ 2 η distribution where η is the difference between the number of estimated parameters in the two models being compared. Standard procedure can then be followed to determine the p-value (p = Pr χ 2 η > t ) for this test and the corresponding significance, where Φ −1 is the inverse cumulative distribution function of the standard normal distribution. Significant results may be tested by either a small p-value or large values of Z. For calculating exclusions we use where L 1 is the likelihood of the best-fit experimental observation. Values of parameters entering θ that result in Z > 2 using this test statistic are then excluded at a 2σ C.L.

VI. RESULTS
We begin with the likelihood maximization test where g x is estimated by fixing the values of both M Z and r. The results are shown in Fig. 4. The 1σ C.L. regions are shaded in the color of the corresponding PDF curves. We repeat the test for two separate values of M Z , 20 MeV [ Fig.4(a)] and 100 MeV [ Fig.4(b)] that represent two choices on either ends of the unconstrained range of values that are preferred by (g − 2) µ . They also show the range of g x favored by the (g − 2) µ anomaly for each value of M Z and the region disfavored by experiments like CCFR [91], Borexino [92], and BaBar [93]. The r = 1 lines correspond to the non-SUSY model (which in this case is the same as the case for degenerate sleptons). The effect of nonvanishing contribution from SUSY sources is captured entirely by choosing other values of r. We show the results for two representative values of r besides the r = 1 result, r = 0.01 and r = 100. These choices capture a wide range of the most common choices for the slepton mass hierarchies.
For all choices of the fixed parameters, we find that the  For perspective, we also show the range of gx that can explain (g − 2)µ (cyan), as well as those values for it that are constrained by experiments like CCFR [91], Borexino [92] and BaBar [93] (gray). Only the unbinned total event counts from the detectors are used for this test. The SM is within 1σ of the best fit value of gx.
range of g x favored by (g − 2) µ is entirely within 1σ of the MLE while the SM is beyond the 1σ C.L. Also, larger values of the stau to smuon mass ratio (r) are slightly preferred over their degeneracy (which is the same as the absence of SUSY in this case) or smaller values. They also have narrower uncertainties, which denotes a better fit to data, and a better overlap with the (g − 2) µ range. All the 1σ intervals of the fit estimates for g x and the deviation of the corresponding MLE from the SM are shown in Table II.
The g x dependence of the chi-square deviation as described in Sec. II D for fixed values of M Z and r = 1 are shown in Fig. 5. The chi-square deviation is defined using total counts only as the full energy-time binned CsI(2020) data is still not in the public domain. The insets in the figures show the variation of the function around its global minima. As in the previous test, we find that the global minima coincides with the (g − 2) µ preferred parameter range for the choices of 20 MeV M Z 100 MeV. The global minima of the chi-square profile is at 3.2 × 10 −4 and 1.3 × 10 −3 for M Z = 20MeV and 100MeV, respectively, which are within the range of couplings allowed from the (g − 2) µ observations. This implies that the band of gauge boson masses from 20−100 MeV better fits the observations of the COHERENT Collaboration than the SM while also explaining (g − 2) µ in tandem. There is an additional minima for larger values of g x that becomes the global minima for M Z 200 MeV but these choices are already ruled out by observations from the experiments mentioned earlier. In contrast to the likelihood profiles, the χ 2 test does not exclude the SM at the 1σ C.L. This distinction is most important as the likelihood test uses both energy and timing distributions from both detectors while the χ 2 test uses only the total counts. However the latter also includes the latest CsI(2020) data while the former does not. Whether the full updated data in a likelihood analysis still has the preference for BSM physics will be clear once the Collaboration makes the data public.
We also show the results of the likelihood maximization test where both g x and M Z are estimated parameters and the ratio r is fixed. The results of this analysis are summarized in Fig.6 where the three different panels demonstrate the effect of three different choices of r. The color code denotes the variation of the PDF. Since the MultiNest algorithm rejects parameter points over iterations by importance sampling, regions of lower likelihood also have a lower density of final reported likelihood points. The red shaded region in these figures denote the parameter space preferred by the COHERENT data within 1σ. We also show the (g − 2) µ preferred and already constrained parameter space in these figures.
There are a few features of these results that are common across all the choices of r. The SM (g x = 0) is always within 2σ of the MLE, however there is a slight preference for nonzero L µ − L τ gauge coupling. Although the MLE is always in regions of the parameter space that have already been excluded by experiments like CCFR [91], Borexino [92] and BaBar [93], unconstrained regions of the parameter space where the (g −2) µ may be explained are within 1σ of it. Changing the value of r affects the estimated most likely value of g x . Larger values of r pre- We also show the parameter space that can explain (g − 2)µ (blue), as well as the region constrained by experiments like CCFR [91], Borexino [92] and BaBar [93] (gray).
fer smaller values of g X and vice versa. Larger values of r also bring the preferred values of g x into better overlap with those favored by (g − 2) µ . At the same time, larger values of this ratio are also associated with a slightly larger likelihood. However, the logarithmic dependence of on r ensures that these dependencies are quite mild as we go for even larger values than 100.
The exclusions from these analyses using the method outlined in Sec. V are shown in Fig. 7. The exclusions from the CENNS10 analysis are far stronger than those from CsI(2017). This is because it observed an excess of events above the SM while the (L µ − L τ ) Z serves to reduce the count for the relevant parameter space. The CsI(2017) data alone has a ∼ 2σ preference for a Z coupling exclusively to second generation leptons [34], leading to a smaller CEνNS count than the SM prediction. Although this preference has a marked reduction in the combined data, it does survive. As a result, the exclu-sions from the combined data closely follow those from CsI(2017) alone. These exclusions have been shown keeping the ratio r = 1 fixed. Since larger ratios show a preference for smaller values of g x , they also correspond to stronger exclusions. However, as mentioned earlier, the dependence of the kinetic mixing on this ratio is logarithmic. Hence, the change in the exclusions get smaller and smaller for larger values of r. The 95% confidence level exclusion from the chi-square profile on the g x -M Z plane is also shown in Fig. 7. These constraints are stronger than the existing ones from the CEνNS detection by the Collaboration, however they are not the strongest in the relevant parameter space.
In addition we also show the projected exclusion from the CCM experiment. Since we have no information for the backgrounds for this experiment yet, we assume that the CEνNS count itself has a Poisson distribution in each bin. We also use a detector threshold of 25 keV and shown. The exclusions shown are for r = 1, which also corresponds to the non-SUSY case. Larger values of r would lead to stronger constraints. The exclusions on the same parameter space from the analysis of the chi-square deviation using the CENNS10 and updated CsI(2020) total counts are also shown.
above that we use 100% efficiency and a quenching factor of 1PE/keV for easy scalability once the information is available. We have used signal predictions using a projected live-period of ∼ 5000hrs/year for three years which correspond to ∼ 3 × 10 22 protons on target. We modify the likelihood function in Eq. 23 suitably by setting the background counts to zero. We also have no information on the associated systematic uncertainties for it, and hence we do not use the dummy variables for marginalization over errors. The only possible error that could be included is the ∼ 2% theoretical error on the CEνNS prediction. However these are too small to cause any significant change to our results and hence are disregarded. Assuming that the CCM experiment observes a CEνNS count equal to the SM prediction, it could potentially exclude hitherto unconstrained portions of the L µ − L τ parameter space. More importantly, it may be able to exclude a portion of the parameter space favored by the (g−2) µ observations. Similar sensitivities in constraining neutrino NSI from the CCM experiment have also been observed by the authors of Ref. [94].

VII. CONCLUSION
We undertook a comprehensive statistical analysis of the combined CsI and CENNS10 data recently made public by the COHERENT Collaboration. The full energy and timing data from the CsI detector as well as the energy, timing and PSD data from the CENNS10 detector was used to estimate the best-fit parameters for the L µ − L τ model and their deviation from the SM prediction. We have also studied the effect that the presence of SUSY may have on these results.
We find a slight preference for reductive BSM effect from the chi-square deviation analysis involving total counts of the updated CsI and CENNS10 data from the COHERENT Collaboration. The global minima of the chi-square profile for 20MeV M Z 100MeV are for values of coupling g x that can explain (g − 2) µ for the respective Z masses and are also unconstrained. An additional minima for larger values of g x is distinct especially for a heavier Z . This additional minima becomes the global minima for M Z 200 MeV, however the choice of this parameter space is constrained by the observations of the CCFR and BaBar Collaborations. A ∼ 1.3σ preference for reductive BSM effect as well as an overlap with the (g − 2) µ favored parameter space within 1σ was also observed from the likelihood profile of the CsI(2017)+ CENNS10 energy, timing and PSD distribution. An important distinction between the two different analyses is that the likelihood profile has no overlap with the SM at the 1σ C.L. In contrast, the χ 2 test using only the total counts from CsI(2020)+CENNS10 data has an overlap with the SM within 1σ despite the global minima mildly preferring nonzero new physics contribution.
The exclusions from the CENNS10 data that observed an excess of counts over the SM CEνNS prediction are much stronger than the earlier CsI exclusions on the L µ − L τ parameter space. However, the surviving preference for reductive BSM effects results in the exclusions from the combined data to closely adhere to the CsI ones. The exclusions from the combined CENNS10 and CsI(2020) total counts using a χ 2 -test are stronger than the CsI(2017) + CENNS10 exclusions, but are still weaker than those coming from the CENNS10 data alone. Still, neither the CENNS10 exclusions by themselves or those from the combined data are the strongest exclusions as far as the L µ − L τ model is concerned. We also show the projected exclusions from the CCM experiment. Assuming it observes the SM prediction for CEνNS , it might be able to exclude yet unconstrained regions of the parameter space favored by the (g − 2) µ .
Some of the other upcoming CEνNS experiments also look promising in this regard. For instance, the planned 610 kg extension (CENNS610) of the CENNS10 LAr detector at the SNS would be able to significantly improve statistics on the existing data. The experiments at ESS and especially ESS10 [95] with its strikingly low threshold energy are expected to probe large parts of the model parameter space including those that can explain the present (g − 2) µ measurements. Other experiments like the NA64µ [96,97] muon beam experiment, the NA62 [98] Kaon factory and liquid Xenon direct Dark Matter detectors like the LZ [99], XENONnT [100] and the DAR-WIN observatory [101] also have promising capabilities to probe regions of the U (1) Lµ−Lτ parameter space pertinent to the (g − 2) µ . For a recent study of projected constraints from these experiments on the L µ − L τ pa-rameter space see, for example, Ref. [102]. Observation of neutrino tridents at the DUNE near detector is also projected to be able to improve upon the CCFR observations and probe the (g − 2) µ favored parameter space of the L µ − L τ model [103,104].
The formalism for the SUSY scenario was kept the same as that in Ref. [80]. Our analysis finds that the range of g x favoured by (g −2) µ,e in this scenario overlaps with the estimates of g x from the likelihood maximization within 1σ. We also find that larger values of the ratio r have better overlap with this range of g x and is slightly more favored by the combined data.
This exercise also points towards the fact that the combined COHERENT data still prefers BSM effects that serve to reduce the CEνNS count predicted from the SM, although to a lesser degree than with the CsI data alone. The CsI(2017) data preferred reductive BSM effects via light mediators at ∼ 2σ [34]. The preference for similar effects in the combined data from CsI(2017) + CENNS10 has been reduced to ∼ 1.3 − 1.4σ.