Nonstandard neutrino interactions at COHERENT, DUNE, T2HK and LHC

We study how nonstandard neutrino interactions (NSI) may be probed by a combination of coherent elastic neutrino-nucleus scattering, neutrino oscillation and collider data, from COHERENT, DUNE, T2HK and the high-luminosity (HL) LHC. We focus on NSI induced by a new flavored gauge boson $Z'$ in a generic anomaly-free ultraviolet-complete model. For $Z'$ masses above 10 GeV, the HL-LHC has the best sensitivity regardless of the flavor structure of the model. For masses between 0.01 GeV-10 GeV, current LHCb data and future COHERENT data have the best sensitivity unless the $Z'$ couplings to the first and second generation leptons are suppressed, in which case DUNE and T2HK have the best sensitivity. For $Z'$ masses between about 5 MeV-20 MeV, DUNE and T2HK have the best sensitivity. We also show how joint analyses of COHERENT and LHC data may constrain such models.


Introduction
Neutrino oscillations have been confirmed by many neutrino experiments using solar, atmospheric, reactor, and accelerator neutrinos in the last two decades. Since the explanation of neutrino oscillations requires nonvanishing neutrino masses, the observation of neutrino oscillation provides clear evidence of new physics beyond the Standard Model (SM) [1]. The next generation precision neutrino oscillation experiments, DUNE and T2HK, will have the sensitivity to probe new physics beyond the standard three neutrino paradigm. A model-independent way of studying new physics in neutrino oscillations was first formulated in Ref. [2], and is now generalized in the framework of an effective field theory for nonstandard interactions (NSI); for reviews see Ref. [3][4][5]. In this framework, NSI not only affect neutrino propagation in matter via neutral current interactions, but also affect neutrino production and detection via charged current interactions. Since modelindependent bounds on the charged current NSI involving charged leptons are generally an order of magnitude stronger than the neutral current NSI [6], we neglect charged current NSI in this work, and focus on neutral current NSI.
In general, neutral current NSI can be described by dimension-six four-fermion operators of the form [2,7], where α, β label the lepton flavors (e, µ, τ ), f denotes the fermion fields (u, d, e), and C indicates the chirality (L, R). Here, 2) with f L αβ , f R αβ being dimensionless parameters that quantify the strength of the new interactions in units of the Fermi constant, G F ≡ ( √ 2v 2 EW ) −1 , with v EW = 246 GeV, the electroweak scale. These contact interactions arise as a result of integrating out a vector mediator significantly heavier than the typical momentum transfer of the processes. As such, the dimensionless coupling parameters are naturally of the order of ∼ g 2 v 2 EW /M 2 , where M and g are the mediator's mass and coupling. Similar to the standard matter effect [2,8], neutral current NSI affect neutrino propagation in matter via coherent forward scattering, in which the momentum transfer is negligibly small compared with other relevant scales involved. Therefore, the adoption of effective four-fermion interactions in Eq. (1.1) is well justified regardless of the mass of the mediator that induces NSI. Also, for neutrinos propagating in unpolarized matter at rest, only the vector combination contributes to the matter potential.
Coherent elastic neutrino-nucleus scattering (CEνNS), which was first observed by the CO-HERENT experiment in 2017 in a cesium iodide (CsI) scintillation detector [9], provides another sensitive probe of new vector neutral current interactions. CEνNS occurs when the momentum transfer Q during neutrino scattering off a nucleus is smaller than inverse of the nuclear radius R. In the process, the scattering amplitudes of the nucleons inside a nucleus are in phase and add coherently, which leads to a large enhancement of the cross section. In the SM, CEνNS is induced via the exchange of a Z boson [10]. Hence, CEνNS is also sensitive to NSI induced by a new neutral vector boson [11,12]. To probe NSI at higher energies, the formulation of Eq. (1.1) may no longer be valid. First, the momentum-dependent propagator of the mediator should be used if its mass M is not much larger than the typical momentum transfer Q to properly model the energy dependence. Second, SM gauge-invariant operators must be adopted when the momentum transfer or the mediator mass is at or above the electroweak scale. Thus, an underlying model that generates NSI is often required. In fact, most of the new physics scenarios associated with the lepton sector at high energies yield NSI [5,6,13].
In this paper we focus on a simple model in which the NSI is induced by a gauge boson Z associated with a new U (1) symmetry. Assuming the presence of three right-handed neutrinos, the most general anomaly-free U (1) model can be generated by with the quark charges Q 1,2,3 and lepton charges Q e,µ,τ satisfying the constraint [14] 3( We further require Q 1 = Q 2 = Q 3 = Q q to avoid large flavor changing neutral currents in the quark sector. The Lagrangian can be written as where the current 1 1 We have decoupled νR assuming they are heavy and inaccessible. with g being the U (1) coupling constant. Since neutrino oscillations are not affected by flavor universal NSI, here we only consider nonuniversal flavor-conserving NSI. Also, because scenarios involving L e are heavily constrained in the low-mass region by electron beam-dump experiments [15][16][17][18][19][20], we set Q e = 0 and only consider the less constrained eletrophobic NSI. For the sake of illustration, we take the following three cases for our benchmark studies [21]: Note that in all three cases the new gauge boson couples to quarks universally. The partial decay width to a pair of fermions is given by where N q = 3, N l = 1, and N ν = 1/2. The branching fractions can then be calculated assuming that the total decay width of the Z is the sum over the SM fermion final states given in Fig. 1.
It is important to note that a SM gauge-invariant formulation of NSI often leads to simultaneous couplings to charged leptons due to the symmetry nature of the gauge doublet 2 (ν, ). This opens up new avenues to search for the new physics associated with NSI, and it also results in stringent constraints on NSI owing to the correlation with the charged leptons. As such, the new gauge boson, if heavy, can be most conveniently searched for at high-energy colliders, especially at the LHC in the di-lepton final state, where X denoted everything in an inclusive search. For our benchmark choices, we have = µ for Cases A and B, and = τ for Case C. We note that in Cases A and B, where muon number L µ is involved, one also can make use of e + e − /pp → 4µ decays at the B-factories and LHC to search for a relatively low mass gauge boson. We do not consider Z bosons lighter than 5 MeV to avoid affecting big bang nucleosynthesis. Once a signal for new physics is observed, it is ultimately important to seek other complementary signals to establish a consistent picture of the underlying physics. In this paper we set out to consider correlated signatures between CEνNS and collider searches.
The rest of the paper is organized as follows. In Section 2, we discuss the current and future sensitivities to NSI from neutrino oscillation experiments. In Section 3, we analyze the current and projected constraints on NSI from the COHERENT experiment. In Section 4, we study constraints on the model from LHC searches. Correlated studies are presented in Section 5. We summarize our results in Section 6.

NSI in neutrino oscillation experiments
The Hamiltonian for neutrino propagation in the presence of neutral current NSI is   Table 1: 2σ allowed ranges for the diagonal NSI parameters from the global analysis of current neutrino oscillation data [24], and from a simulation of DUNE and T2HK.
and V is the potential from interactions of neutrinos in matter, which can be expressed using the NSI operators in Eq. (1.1) as Here, V CC ≡ √ 2G F N e , is the standard matter potential, and the effective NSI parameters are with N q,e the number density of fermions q = u, d and e.
Since neutrino propagation in matter is affected by coherent forward scattering, in which the momentum transfer is zero, the effective Lagrangian from Eq. (1.5) that is relevant for NSI can be written as We can then use the bounds on the NSI parameters from neutrino oscillation experiments to constrain the parameter spaces in the Z models. For Case A (C), the model predicts that only µµ ( τ τ ) is nonzero. For Case B, since µµ is equal to τ τ , and neutrino oscillation probabilities are not affected by a subtraction of a diagonal contribution from the full Hamiltonian, we can obtain constraints on Case B from bounds on NSI with only ee being nonzero. We adopt the 2σ bounds on u αα from the global analysis of current oscillation data [24] as compiled in Table 1. Note that neutrino oscillation data constrain differences between two diagonal 's, not individual diagonal 's. To obtain bounds on a single , we set one of the two 's to be zero. We bound u µµ by choosing the smaller of the values obtained by setting u ee = 0 in u ee − u µµ and u τ τ = 0 in u τ τ − u µµ . We apply them to constrain the theory parameter space in the (M Z , g ) plane and the exclusion regions are shown as the purple areas in Fig. 2. Note that the bounds from the global analysis are obtained under the assumption that all NSI parameters are nonzero and then projected to one NSI parameter. Since degeneracies among NSI parameters . The red shaded areas correspond to the 2σ exclusion regions by using the energy spectrum from the COHERENT CsI detector [9]. The red dashed lines show the expected 2σ limit from COHERENT with a 750 kg LAr detector [25] and a 4-year exposure using both energy and time information. The purple areas correspond to the 2σ bounds from a global fit to neutrino oscillation data [24]. The dashed purple lines show the expected 2σ exclusion limit from DUNE and T2HK combined. Regions above the brown curves are excluded by CMS [26] and BaBar [27] at 2σ and 90% CL, respectively, using pp/e + e − → µ + µ − Z searches. The brown dashed curves are the 2σ expected sensitivities from HL-LHC, with an integrated luminosity of 3000 fb −1 , in the µ + µ − Z channel, and the blue solid (dashed) curves correspond to the expected 2σ (5σ) limit using di-muon searches for Cases A and B, and di-tau searches for Case C. In the upper panels, the blue shaded regions are excluded at 90% CL by the LHCb dark photon searches [28] and at 2σ by the ATLAS di-muon searches [29] with 139 fb −1 . In the lower panel, the blue area is excluded at 2σ by the ATLAS di-tau searches [30] with 36.1 fb −1 . The 2σ limit from CCFR [31,32] is given by the orange curves. The 2σ allowed regions that explain the discrepancy in the anomalous magnetic moment of the muon (∆a µ = (29 ± 9) × 10 −10 [33]) are indicated by the black band. The black stars mark the benchmark points we consider in Section 5.
can significantly weaken the constraints on an individual NSI parameter [34], the current bounds from the global analysis of oscillation data should be considered to be conservative.
We also consider the sensitivity of the next generation long-baseline neutrino oscillation experiments, DUNE [35] and T2HK [36]. We follow the procedure of Ref. [37], and simulate the DUNE and T2HK data assuming the normal neutrino mass hierarchy, the neutrino CP phase δ = 0, and αα = 0. We scan over both the mass hierarchies, the neutrino oscillation parameters and take only one diagonal αα to be nonzero at a time. The 2σ allowed ranges for the diagonal NSI parameters are provided in the last column of Table 1. The expected sensitivities in the (M Z , g ) parameter space are shown as the purple dashed lines in Fig. 2. As expected, it simply scales linearly with g /M Z . The reaches for the three cases are roughly similar. For instance, at M Z ∼ 10 GeV, the sensitivity for the couplings can reach g ∼ 0.008 (0.02) [0.008] for Case A (B) [C]. We see that future bounds on NSI will be improved by a factor of a few compared to current bounds, and the current constraints on the parameter space in Case C for M Z 200 GeV only come from neutrino oscillation data.

CEνNS
CEνNS has recently been measured by the COHERENT experiment, which detects neutrinos from the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory. Neutrinos at the SNS [38] consist of a prompt component of monoenergetic ν µ from the stopped pion decays, π + → µ + + ν µ , and two delayed components ofν µ and ν e from the subsequent muon decays, µ + → e + +ν µ + ν e . The fluxes of the three neutrino flavors (ν µ ,ν µ and ν e ) are well known and given by where the normalization factor N = rT N POT 4πL 2 , with r = 0.08 the number of neutrinos per flavor produced per proton collision, N POT = 2.1 × 10 23 the total number of protons on target per year, T the number of years of data collection, and L the distance between the source and the detector [9]. The ν µ component has energy (m 2 π − m 2 µ )/(2m π ) ≈ 30 MeV, and the energies of the ν e andν µ have a kinematic upper bound, m µ /2 ≈ 50 MeV. The expected number of events with recoil energy in the energy range [E r , E r + ∆E r ] and arrival time in the time interval [t, t + ∆t] is given by where m det is the detector mass, M is the molar mass of the target nucleus, N A = 6.022 × 10 23 mol −1 , ρ α (t) is the arrival time Probability Density Function (PDF) provided in the CO-HERENT data release [39], and α = ν µ ,ν µ , ν e . We assume that the presence of new neutral current interactions do not modify the arrival time PDF.
Neglecting radiative corrections, the differential cross section for a given neutrino flavor ν α scattering off a nucleus is given by where F (Q 2 ) refers to the nuclear form factor taken from Ref. [40]. In the presence of NSI, the effective charge can be written as where Z (N ) is the number of protons (neutrons) in the nucleus, g V p = 1 2 − 2 sin 2 θ W and g V n = − 1 2 are the SM weak couplings, and θ W is the weak mixing angle. The NSI parameters for coupling to up and down quarks can be written as .

(3.5)
For the CsI detector, the total cross section is a sum of the contributions of 133 Cs and 127 I, i.e., To compare with COHERENT data, we convert the nuclear recoil energy to the number of photoelectrons (n PE ) using the relation [9], Note that we do not use the new quenching factor reported in Ref. [41] as it is still under investigation by the COHERENT collaboration [42]. We employ the acceptance function [39], where k 1 = 0.6655, k 2 = 0.4942, x 0 = 10.8507 and θ(x) is the Heaviside step function. Because the number of events is small and experimental uncertainties large, we use the energy spectrum (but not the timing information) measured by the CsI detector to evaluate the statistical significance of a nonstandard scenario. We define 9) where N i meas and N i th is the number of measured and predicted events per energy bin, respectively. The statistical uncertainty per energy bin is σ i where B SS and B on are the estimated steady-state and beam-on backgrounds, respectively. B SS is determined by the anti-coincident (AC) data, and B on mainly consists of prompt neutrons. Both the spectral and temporal distributions of the backgrounds are provided by the COHERENT collaboration [39]. For the signal normalization uncertainty, we follow the original COHERENT analysis and choose σ γ = 0.28, which includes the neutrino flux uncertainty (10%), form factor uncertainty (5%), signal acceptance uncertainty (5%), and quenching factor uncertainty (25%). For the beam-on background uncertainty, we fix σ β = 0.25 [9]. We scan over values of the coupling g and the mediator mass M Z . The 2σ exclusion regions in the (M Z , g ) plane are shown as the red regions in Fig. 2 for Cases A and B. For M Z 50 MeV, the current constraint from COHERENT CsI is comparable to the expected sensitivity of DUNE+T2HK for Case B, and is weaker by about a factor of two for Case A. For very small M Z DUNE+T2HK has greater sensitivty than the current COHERENT bounds for both Cases A and B. Note that COHERENT data does not place bounds on Case C because the SNS beam does not have ν τ andν τ .
The COHERENT collaboration has an extensive upgrade plan [25], part of which is a 750 kg LAr detector located at L = 29 m from the source. We assume a 4-year exposure with the same neutrino production rate as the current setup, which corresponds to 8.4×10 23 protons-on-target (POT) in total. Since both the spectral and temporal distributions of the recoil energy events depend on the flavor structure, we perform a two dimensional analysis that utilizes both the spectral and temporal information. To estimate the projected sensitivities at the LAr detector, we adopt the likelihood function from Ref. [43], i.e., where λ(t, E r ) = (1+γ)N th (t, E r , )+βN obs,bg (t, E r ). We calculate the number of events expected in the SM for each bin within the range 0 < t < 6 µs and 20 keV < E r < 100 keV, with bin sizes of 0.5 µs and 2 keV, respectively. We assume that the steady-state background is uniform in energy and is 1/4 of the SM expectation. We also assume the systematic uncertainty σ γ to be 17.5%, which corresponds to a reduced quenching factor uncertainty of 12.5% for LAr. A more precise treatment would include energy-dependent form factor uncertainties [44]. The projected sensitivities are shown by the purple dashed line in Fig. 2. A factor of three improvement is expected in the sensitivity to the coupling, compared to the current CsI results. We see that future CEνNS experiments will set stronger bounds than next generation neutrino oscillation experiments for most Z masses in Cases A and B, and will provide the strongest constraints for 20 (10) MeV M Z 1 GeV in Case A (B).

Collider searches for NSI
As emphasized in the introduction, a SM gauge-invariant formulation of NSI often results in simultaneous couplings to charged leptons. This opens up new avenues to search for the new physics associated with NSI, in particular at colliders. We explore the sensitivity reach at the LHC for NSI via a di-lepton final state from the Drell-Yan (DY) production of a Z , with = µ, τ and X denotes other inclusive states (like a jet) when kinematically favorable for the signal identification. This is a particularly sensitive signal M Z > M Z . We also include a four-lepton final state, This channel is more suitable for a low mass Z as we will see below.
We use the Monte Carlo event generator MadGraph5 aMC@NLO [45] to generate signal and background samples with the NN23LO1 PDF set [46]. The NSI Lagrangian is implemented in the FeynRules 2.0 [47] framework. Pythia 8.1 [48,49] is used for parton showering and hadronization. Matching is performed with the MLM prescription [50]. The generated events are passed into Delphes 3.4.1 [51] for fast detector simulation.

Cases A and B: µ final states
In Case A, the new gauge boson couples to quarks universally, and only to second generation leptons. While in Case B, the new gauge boson couples equally to second and third generations leptons. We first apply the existing LHC bound on searches for the di-muon final state to both cases, given that muons are much easier to identify than taus at the LHC. ATLAS [29] has performed a search for di-lepton resonances in the 250 GeV M Z 6 TeV mass range setting a 2σ upper limit on the fiducial cross section times branching ratio with 139 fb −1 at √ s = 13 TeV. The fiducial region is defined by the acceptance cuts, To extract limits on g , we calculate σ(pp → Z + X) · B(Z → µ + µ − ) in the fiducial region at leading order (LO). The expected signal yields are rescaled to next-to-leading order (NLO) accuracy using a K-factor of 1.3 [52]. From the auxiliary figure 2c of Ref. [29], the upper limits at 2σ on the fiducial cross section from ATLAS are translated into the bounds on our model parameters, shown as the blue shaded regions in the upper panels of Fig. 2. This search excludes g 1.6 (2.4) × 10 −3 for M Z ≈ 250 GeV in Case A (B). Searches for dark photons decaying to di-leptons can shed light on new vector bosons, especially relatively light ones. In Cases A and B, we recast prompt-like dark photon searches at LHCb [28] to obtain constraints in the mass range 200 MeV to 70 GeV based on the framework developed in Ref. [53]. This is the most sensitive probe currently in this mass window except near the resonances like J/ψ, Υ and approaching the Z-pole. The corresponding upper limits on the coupling at 90% CL are shown by the blue shaded regions in Fig. 2.
Having discussed the bounds from the di-muon final state, we turn to the four-muon final state. Both the BaBar and CMS have performed searches for the decay, γ * /Z * → µ + µ − Z → 4µ. The BaBar searches [27] set a 90% CL upper limit on the new gauge coupling based on a L µ − L τ model corresponding to Q q = Q e = 0, Q µ = −Q τ = 1 in our parameterization. The CMS searches [26] set a 2σ upper limit on g by assuming the branching ratio B(Z → µ + µ − ) = 1/3 and Q µ = 1. By rescaling the observed bounds according to the branching fractions and production cross section, we extract bounds for our scenarios. The brown curves show the BaBar and CMS bounds in the upper panels of Fig. 2. We see that the current bound from the LHCb dark photon search is dominant in the medium mass range and disfavors g 10 −4 for M Z ≈ 200 MeV.
We further estimate the sensitivity reach via the di-muon channel Z → µ + µ − for 10 M Z 6000 GeV at the high-luminosity LHC (HL-LHC) with the full 3000 fb −1 integrated luminosity. The signal is from the DY process as in Eq. (4.1). We select events that contain at least two opposite-sign muons. The leading (subleading) muon is required to have p T > 22 (10) GeV. All muons are required to have |η| < 2.4. Finally, in calculating the sensitivity, we apply a mass window cut 0.97 M Z < M ( + − ) < 1.03 M Z below 3 TeV, and use a 3 − 6 TeV mass window to ensure enough background events in the high mass region, to optimize the signal observability. The dominant background is from the SM DY process. We also include smaller background contributions from tt, tW , W W and ZZ. We generate the signal and DY background with up to two additional jets in the phase space M µµ < 60 GeV. This is so that for a lighter Z , the additional jets help to kick the leptons to a high momentum for more efficient triggering. For M µµ > 60 GeV, we generate the signal and DY background at LO and apply the combined QCD and electroweak corrections to the invariant mass distributions according to Ref. [54]. tt and tW backgrounds are generated at LO and normalized to NNLO + NNLL by a K-factor of 1.84 [55] and 1.35 [56] respectively. The W W, W Z, and ZZ backgrounds are normalized to NNLO QCD by a K-factor of 1.98 [57], 2.07 [58], and 1.74 [59] respectively. The local significance is defined as 4) where N S (N B ) is the expected number of signal (SM background) events. The blue solid (dashed) curves in the upper panels of Fig. 2 show the 2σ (5σ) sensitivities. The sensitivity is significantly improved in a broad mass range.

Case C: τ final states
For Case C, the signal channel at the LHC is pp → Z + X with Z decaying to a tau pair. For a high-mass mediator decaying to di-tau, ATLAS [30] and CMS [60] have set a 2σ upper limit on inclusive σ(pp → Z + X) · B(Z → τ + τ − ) in the 200 GeV M Z 4 TeV (ATLAS) and 500 GeV M Z 3 TeV (CMS) mass ranges with √ s = 13 TeV and 36.1 fb −1 and 2.2 fb −1 , respectively. We only display the ATLAS constraint on g in the lower panel of Fig. 2.
We also estimate the sensitivity reach for 20 GeV M Z 6000 GeV at the HL-LHC with 3000 fb −1 of integrated luminosity. There are mainly four decay modes for di-tau, namely, τ e τ µ (6%), τ e τ h (23%), τ µ τ h (23%), and τ h τ h (42%), where h denotes a hadron. In this analysis, we use the TauDecay package [61] to model the relatively clean leptonic and semi-leptonic decay modes of the taus. The main backgrounds for τ e τ µ are tt, W W , and DY. For the semi-leptonic modes, the main backgrounds are DY and W +jets. To include the QCD multijet background in the semi-leptonic modes, we add 6% and 28% of the sum of the DY and W +jets backgrounds for the τ µ τ h and τ e τ h modes, respectively [60]. The signal and DY background events are generated at LO and scaled by a K-factor of 1.3 [52] for M τ τ > M Z , while for M τ τ < M Z , we generate the signal and DY background with up to two additional jets in the final states. We generate tt, W W , and W +jets background events at LO. To take higher-order corrections into account, the LO cross section of tt is normalized to the NNLO + NNLL cross section by a factor of 1.84 [55]. The LO cross sections of W W and W +jets are normalized to NNLO QCD by a factor of 1.98 [57] and 1.46 [62], respectively. To reduce the background, we implement two different selection rules SR1 and SR2 for M Z below and above the Z-pole. In the τ e τ µ mode, both SR1 and SR2 require: • Only one muon and one oppositely charged electron with p T > 20 GeV and |η| < 2.4, • veto b-tagged jets, where τ 1 and τ 2 are respectively e and µ, and M µ T is the transverse mass of the charged lepton µ and the missing transverse momentum / E T is defined as In addition, SR1 requires where ∆R is the angular distance between τ 1 and τ 2 . ∆R cut is varied with M Z to maximize the local significance S l . For example, we choose ∆R cut = 1.0 (1.6) for M Z = 20 (40) GeV. SR2 further requires • cos ∆φ(τ 1 , τ 2 ) < −0.95 , (4.6) • cos ∆φ(τ 1 , / E T ) + cos ∆φ(τ 2 , / E T ) > −0.1 , (4.7) where the missing energy cut / E cut T is varied with M Z to maximize the local significance S l . We take / E cut T to be 40 (450) GeV for M Z = 500 (2000) GeV. In the τ τ h modes, both SR1 and SR2 require: • Only one charged lepton and at least one opposite-sign tau-tagged jet with p T > 20 GeV and |η| < 2.4, • veto b-tagged jets, • M T < 40 GeV. (4.9) The further requirements of SR1 and SR2 are the same as for the leptonic τ e τ µ mode, with τ 1 and τ 2 the charged lepton and tau-tagged jet, respectively. The blue solid (dashed) curve in the lower panel of Fig. 2 shows the 2σ (5σ) sensitivity for Case C using a combination of the three decay modes (τ e τ µ , τ e τ h , and τ µ τ h ), respectively, with 3000 fb −1 at the HL-LHC.

Correlated signatures at CEνNS and collider experiments
It is of fundamental importance that we observe correlated signals of NSI in different experiments. In this section, we study correlated signatures at future CEνNS and collider experiments. We first simulate spectra in the presence of NSI and then examine the consistency between the two experiments in the hope of identifying a correlated signal. We select the benchmark point, for Cases A and B and explore how a signal observed in one experiment will manifest in another. The point is marked with a star in Fig. 2. The point is chosen so that observable signals can be produced at COHERENT and at the LHC. Since this set of parameters does not produce a signal at DUNE and T2HK, we focus on correlated signatures at COHERENT with an upgraded LAr detector and the high luminosity LHC with L = 3000 fb −1 . Note that the benchmark point is chosen in a currently allowed narrow region near m(Υ(1S)), and that LHCb data impose strong constraints for M Z below and above it. We first study signatures at COHERENT with an upgraded LAr detector. The recoil energy and temporal distributions of the events are shown in the left and right panel of Fig. 3, respectively. As can be seen from the left panel, the event excess is mainly at low energies. From the right panel, we see that the event excess peaks at around t = 1 µs. This is due to the fact that the prompt component of the COHERENT flux is primarily composed of ν µ , and the NSI coupling to ν µ leads to a modification of the number of events in Cases A and B. To analyze the spectra and to facilitate a joint analysis with simulated LHC data, we define where L( θ) is defined in Eq. (3.10) with θ = {g , M Z }. We then calculate ∆χ 2 = χ 2 − χ 2 min . The 2σ allowed region for Case A and 1σ allowed region for Case B, with data simulated with our benchmark point, are the regions between the red curves in Fig. 4. The 2σ regions for Case B are too large to display.
We now study signatures at the HL-LHC. Since we are interested in the low-mass region, we focus on the clean channel, Z → µ + µ − Z → 4µ. We generate the leading process qq → 4µ at the leading order (LO). Following the CMS analysis [26], we require at least four well-identified and isolated muons to have p T > 5 GeV and to be in the central region of the detector |η| < 2.4, with at least two muons to have p T > 10 GeV and at least one to have p T > 20 GeV. Dimuon candidates formed from an opposite sign muon pair are required to have 4 < M µ + µ − < 120 GeV.  Figure 4: 2σ allowed regions for Case A (left) and 1σ allowed regions for Case B (right) from COHERENT with a large LAr detector (within the red curves) and HL-LHC Z → 4µ decays (within the blue curves). The purple shaded regions (2σ for case A and 1σ for Case B) are from our joint analysis. The magenta shaded regions are the allowed regions after including the LHCb bound as a prior. The stars mark the best fit points from our joint analysis.
The four selected muons are required to have zero net charge and 80 < M 4µ < 100 GeV. The NNLO/LO K-factor is chosen to be 1.29 [26]. By following the CMS procedure in Ref. [26], we are able to reconstruct M Z , whose distributions are shown in the left panel of Fig. 5. Unfortunately for Z s of GeV mass, COHERENT sees an overall suppression in the CEνNS event rate, but no spectral distortion, thereby precluding it from determining M Z . So a di-muon invariant mass cut cannot be applied and the look-elsewhere effect must be taken into account. Instead, we employ the M 4µ distributions (shown in the right panel of Fig. 5) to evaluate the precision with which the Z parameters can be determined. We divide the range of M 4µ (80 GeV, 100 GeV) equally into 10 bins and perform a χ 2 analysis with where N S,i (N B,i ) is the expected number of signal (background) events in the i th bin. The background systematic uncertainty σ B is chosen to be 5%. The parameters favored at 2σ for Case A and at 1σ for Case B lie between the blue curves in Fig. 4; Case B has no lower blue curve because the SM is allowed at 1σ. (The brown dashed curves in Fig. 2 for the 2σ sensitivity to the 4µ channel are produced by requiring the di-muon invariant mass M µ + µ − to be within 2% of M Z , and defining the local significance as N S / N B + σ 2 B N 2 B .) We perform a joint analysis of future COHERENT and HL-LHC data by combining the two χ 2 in Eqs. (5.1) and (5.2). The resulting 2σ allowed regions for Case A and 1σ allowed regions for Case B are shaded in purple in Fig. 4. Consider Case A. The fact that the allowed regions from COHERENT and LHC have different slopes enables a combination of their datasets to limit M Z to be below about 60 GeV. However, a precise determination of M Z is not achieved even by combining the datasets. For Case B, both COHERENT and HL-LHC only provide upper bounds on g at 2σ. COHERENT dominates the sensitivity and the HL-LHC does not lead to a clear signal observation in the parameter region considered.
We now impose the stringent bounds from LHCb. To include the LHCb constraint, for each value of M Z we add χ 2 LHCb = 2.71(g /g bound ) 2 to our joint χ 2 , where g bound is the 90% CL exclusion limit from LHCb at that value of M Z ; note that the LHCb dark photon search [28] is performed independently at each mass, so that only one parameter, g , is varied in the analysis. On including the LHCb constraint, the allowed regions shrink significantly; see the magenta shaded regions in Fig. 4.

Summary
Next generation neutrino oscillation and CEνNS experiments will reach the sensitivity to discover new physics parameterized in the form of NSI. It is natural to seek complementary probes for NSI. Indeed, most beyond the Standard Model scenarios that generate NSI often result in simultaneous couplings to charged leptons, which opens up new possibilities to search for new physics associated with NSI at colliders.
In this work we studied a simple anomaly-free, ultraviolet-complete, gauged U (1) model that generates lepton flavor universality violating NSI. We considered three scenarios: B − 3L µ , B − 3 2 (L µ + L τ ), and B − 3L τ . The Z decay branching fractions are shown in Fig. 1. Our main results are shown in Fig. 2 and we summarize them as follows. For constraints from current data: 1. In Cases A and B, we mainly use neutrino oscillation, CEνNS, and collider experiments to put constraints on the coupling g in the mass range, 5 MeV< M Z < 6 TeV. We found that neutrino oscillation and CEνNS experiments give the most stringent bounds for masses below the dimuon threshold which is around 200 MeV. Above the dimuon threshold up to 70 GeV, LHCb prompt-like dark photon searches provide the strongest constraints except near the J/ψ, Υ resonances and in the vicinity of the Z-pole. ATLAS dimuon searches give the strongest bounds in the mass range, 250 GeV ≤ M Z ≤ 6 TeV.