$Z^\prime$s in neutrino scattering at DUNE

Novel leptophilic neutral currents can be tested at upcoming neutrino oscillation experiments using two complementary processes, neutrino trident production and neutrino-electron ($\nu-e$) elastic scattering. Considering generic anomaly-free $U(1)$ extensions of the Standard Model, we discuss the characteristics of $\nu-e$ scattering as well as $e^+e^-$ and $\mu^+\mu^-$ trident production at the DUNE near detector in the presence of such BSM scenarios. We then determine the sensitivity of DUNE in constraining the well-known $L_e - L_\mu$ and $L_\mu - L_\tau$ models. We conclude that DUNE will be able to probe these leptophilic models with unprecedented sensitivity, covering unproved explanations of the $(g-2)_\mu$ discrepancy.


I. INTRODUCTION
The discovery of neutrino oscillation is the first laboratory-based proof of physics beyond the Standard Model (BSM) establishing that, in contrast to the predictions of the Standard Model (SM), the neutrino sector has at least three mass eigenstates distinct from the flavor states defined by the charged leptons. However, the mechanism that generates neutrino masses remains unknown and many competing candidate theories exist, ranging from the simplicity of a Dirac mass term protected by a symmetry (see, e.g., [1][2][3]) or the popular seesaw mechanisms [4][5][6][7][8][9][10][11][12] to proposals with a more elaborate spectrum of particles. In general, more elaborate scenarios have additional motivations, including the explanation of lepton mass and flavor hierarchies (see, e.g., [13]), the matter-antimatter asymmetry of the universe [14][15][16], the existence of dark matter [17,18], the scale of neutrino masses [19][20][21], or the anomalous experimental results [22]. Uncovering the nature of new physics in the neutrino sector, and its connection to other BSM concerns, will be a central aim of the experimental and theoretical programs over the next few decades.
Although significant progress has already been made, the neutrino sector remains relatively poorly explored. There are still large uncertainties on the masses and on some of the mixing parameters of the light neutrinos [23], but even beyond the effects of neutrino mass, many SM cross sections are poorly known theoretically and infrequently measured. This is in part due to the typical energy scales of neutrino experiments that necessitate the modeling of the neutrino-nucleus interactions, but also because of the rareness of neutrino scattering (see Ref. [24]). Much effort has gone into measuring crucial cross sections at oscillation experiments [25][26][27] and at the Main Injector Experiment for ν-A (MINERνA) [28], a dedicated cross section experiment. However, given the necessity and potential richness of BSM physics in the neutrino sector, and the wide array of measurements yet to be made, it is conceivable that new physics will also manifest itself as detectable signatures in neutrino scattering. It is crucial to keep an open mind about what future experimental work might find, for instance, in the auxiliary physics program of the near detector (ND) of the next-generation Deep Underground Neutrino Experiment (DUNE) [29].
Novel interactions in the neutrino sector have been proposed for a variety of reasons, including as a potentially observable effect in the neutrino oscillation probabilities (see, e.g., [30]), as a way of ameliorating tension introduced by sterile neutrinos in the early universe [31][32][33][34][35][36][37][38], and as a possible explanation of anomalous results at a short baseline [39][40][41]. Models that introduce new interactions between neutrinos and matter have been discussed in simplified settings [42][43][44], via effective field theory [45,46] and specific UV complete models [47] (see also [48] for a neutrinophilic Z 0 study at the DUNE ND). One class of models restricts the new interactions to leptons. This arises most naturally in settings with a gauged subgroup of the lepton number, with most attention given to the anomalyfree subgroups L α − L β for α; β ∈ fe; μ; τg [49,50]. Such leptophilic interactions must satisfy strong constraints from processes involving charged leptons [51], but in the case of a gauged L μ − L τ symmetry, neutrino processes have been found to be particularly competitive [52].
In this work, we study potential constraints that can be placed on a general set of leptophilic Z 0 models in the two most likely scattering channels for this type of BSM at the near detector of DUNE: ν − e scattering and νll trident scattering. During ten years of running, a 75-t near detector subjected to the intense neutrino beam at the Long-Baseline Neutrino Facility (LBNF) will provide tens of thousands of ν − e scattering events. The cross section for this process is theoretically well understood and can therefore be a sensitive probe of BSM physics. Additionally, this process has received special interest due to its potential in reducing systematic uncertainties in the neutrino flux [53,54], an undertaking that can be affected by new physics. Despite not being a purely leptonic process, neutrino trident production can also be measured with reasonable precision at DUNE, where hundreds of coherent and diffractive trident events are expected at the ND [55]. We study the neutral current channels with dielectron or dimuon final states, pointing out how the new physics contribution impacts the nontrivial kinematics of these processes. The main advantage in such measurements lies on the flavor structure of dimuon tridents, which can be used to constrain otherwise difficult to test models, such as the one where a new force is associated with the L μ − L τ gauge symmetry [52].
Although these processes can place stringent bounds on many classes of mediators, many scenarios are already heavily constrained through other experimental work. A recent study of several different Uð1Þ X models using ν − e scattering was presented in Ref. [56], where data from past ν − e experiments CHARM-II, GEMMA, and TEXONO have been used to put bounds on the couplings and masses of general Z 0 s. Novel charged particles are typically constrained to be very massive, leading to little enhancement of the charged current neutrino scattering rates. In particular, charged scalars have been considered in νll trident scattering in Ref. [57], where it is found that trident measurements can provide competitive bounds on charged scalars, albeit only in simplified theoretical settings. The requirement of doubly charged scalars or the connection to neutrino masses introduced by the typical UV completions of such models dilutes the relevance of the trident bounds. Neutral scalars are viable but also present challenging UV completions. Novel Z 0 interactions in νll trident scattering with dimuon final states have been studied in Ref. [52], where it was shown to be a promising channel to probe a L μ − L τ gauge symmetry. This model was revisited in Refs. [58,59], where the effects of kinetic mixing and the possibility of a measurement by T2K was alluded to. Finally, neutrino trident scattering with atmospheric neutrinos was shown to be sensitive to this model as well as to simplified scalar models in [60]. It should be noted, however, that as it was shown in Ref. [55] the equivalent photon approximation (EPA) discussed in several recent studies [57,58] for the calculation of the trident cross section leads to intolerably large errors in the predictions for the νll scattering channels in the SM. For this reason, we calculate this process without making this approximation.
The structure of the paper is as follows. In Sec. II, we describe the basic properties of the leptophilic scenarios that we will consider in this work. In Sec. III, we discuss the calculation of the trident and ν − e scattering cross sections in a general model of leptophilic Z 0 . In Sec. IV, we show how DUNE can place bounds on a few popular leptophilic Z 0 models, discussing our assumptions for experimental configurations and backgrounds. We make our concluding remarks in Sec. V.

II. LEPTOPHILIC Z 0 MODELS
Since we are interested in models where the novel neutral currents are present only in the lepton sector, let us consider explicitly a Uð1Þ Z 0 extension of the SM whose Lagrangian is given by where L α (l α ) represents the leptonic SU(2) doublet (singlet) of flavor α ∈ fe; μ; τg, and we included N right-handed neutrinos with charges Q N under the new symmetry for completeness. Thus, we have 7 þ N new parameters to characterize the couplings between the new boson and the lepton sector, one gauge coupling g 0 and 6 þ N charges fQ L α ; Q R α ; Q N g. Below the scale of the electroweak symmetry breaking (EWSB), the relevant interaction terms in the Lagrangian are given by We note that the right-handed singlets could modify the form of the neutrino interaction in Eq. (2) by introducing a right-chiral current. The details of this would depend on the relationship between these chiral states and the flavor-basis neutrino ν α . However, in practice our Lagrangian is fully general, as the polarization effects in the neutrino beam ensure that only the left-handed charge is relevant for lightneutrino scattering experiments.
The Lagrangian in Eq. (1) contains all of the terms necessary for this analysis. However, when it comes to assigning specific charges to the particles, a few wider model-building considerations are worthy of discussion. In the SM, any nonvectorial symmetry would forbid the Yukawas responsible for the charged-lepton mass terms post-ESWB; similarly, possible negative implications for neutrino mass generation are expected. The precise implementation of the neutrino mass mechanism is highly model dependent, but neutrino gauge charges are not compatible with many usual realizations. 1 Furthermore, the novel gauge boson Z 0 will also require a mass generation mechanism, and indeed this could be achieved via the means of symmetry breaking. Although each of these is an important aspect of model building, their resolution can be expected to have little impact on the phenomenology of neutrino scattering, and we will not pursue them here. Anomaly freedom of our new symmetry, however, is a more pertinent concern. It has been shown that an anomalous group can always be made anomaly free via the introduction of exotically charged sets of fermions which can be given arbitrarily large masses [61]. Yet these novel fermionic states necessarily introduce effects at low scales, which in some cases can strongly affect the phenomenology of the model [62]. Therefore, while it seems likely that mass generation can be addressed with the addition of new particles that do not interfere with neutrino scattering phenomenology, anomaly freedom is more pernicious. For this reason we will briefly discuss how anomaly freedom will dictate the types of leptonic symmetries that we consider in the remainder of this work.
a. Anomaly freedom. The most general anomaly-free symmetries compatible with the SM were first deduced in the context of grand unification theories [63,64]. More recently, an atlas of all anomaly-free Uð1Þ extensions of the SM with flavor-dependent charges has been provided by Ref. [65]. Interestingly, the only anomaly-free subgroups of the SM with renormalizable Yukawa sectors are leptophilic: the lepton-family number differences L α − L β ðα; β ¼ e; μ; τÞ [49,50]. The popular B − L symmetry is, in fact, anomalous unless right-handed SM singlets are added with the appropriate charges. This is well motivated by the necessity of neutrino mass generation but remains a hypothesis, as not all models of neutrino mass require novel fermionic content. For the sake of discussion, we follow a similar logic and consider the most general anomalyfree subgroups of the SM accidental leptonic symmetries allowing for an arbitrary number of right-handed fermionic singlets. These would presumably be associated with the neutrino mass generation mechanism, but we impose no specific relations in this regard due to the significant modelbuilding freedom. The anomaly conditions for a leptophilic model with right-handed neutrinos are given below 2 [66]: Gauge-gravity In the absence of new N R particles (Q N ¼ 0) and assuming that Q L α ¼ Q R α , that is, considering vector couplings, we find the three well-known discrete solutions for Eqs. (3): the antisymmetric pairs L α − L β , α; β ¼ fe; μ; τg; α ≠ β. As far as anomalies are concerned, all three pairs are equal, but frequently focus falls on L μ − L τ , which has no coupling to electrons and correspondingly weaker constraints. If we reconsider these conditions with charged right-handed neutrinos, we find a one-dimensional continuous family of potential symmetries that can be consistently gauged. We can parametrize this as 1 If neutrino masses are thought of as coming from a Weinberg operator, it is clear that the leptonic doublet must be uncharged under any unbroken Uð1Þ 0 group. 2 Notice that Uð1Þ 3 Z 0 together with gauge-gravity conditions imply that the number of right-handed states must be at least What we have shown is that linear combinations of the ðL α − L β Þ choice of charges yield an anomaly-free scenario provided N right-handed neutrinos respecting Eq. (5) are added to the theory. We have checked that the "anomalyfree atlas" in [67] contains a subset of these solutions, which are more general.
The above conclusions are based on the assumption of vectorial charge assignments. In the SM, this requirement is a consequence of the origin of mass assuming a chargeless Higgs. However, in nonminimal models this requirement could be relaxed. Even with this extra freedom, not all charge assignments are allowed: for example, a purely chiral Uð1Þ 0 cannot satisfy Eq. (3c) without additional matter charged under the SM gauge group. The axial-vector case, however, does have further solutions: we find that the same onedimensional family of charges is allowed as for the vectorial gauge boson-in this case, the charges apply to the lefthanded fields and the right-handed ones have the opposite charges. In such a model the leptonic mass generation mechanism would necessarily be more complicated than in the SM, but such a possibility is not excluded. UV completions of an axial-vector Z 0 have been presented in [68,69]; however, these generally introduce extra bounds that are expected to be stronger than neutrino scattering bounds (see, e.g., [62,70]). For this reason, we only comment on the consequences of an axial-vector case in our calculations, but do not develop any particular model or constraint.
b. Kinetic mixing. The symmetries of our SM extensions allow for kinetic mixing between the Z 0 and the SM gauge bosons [71][72][73] where F κρ and F 0κρ are the field strength tensors of the hypercharge and the Z 0 boson, respectively. The presence of such coupling introduces a very rich phenomenology and has been explored in great detail in the literature [74]. In this work, we choose to focus on the less constrained possibility of vanishing tree-level kinetic mixing. In this case, kinetic mixing is still radiatively generated due to the presence of particles charged under both the SM and the new Uð1Þ group. As well as the SM particle content, additional particles present in the UV theory may also contribute to kinetic mixing, but we will neglect these contributions in this study as they are highly model dependent. 3 We compute ε between the Z 0 and the SM photon, and find the one-loop result to be finite for any ϱðL α − L β Þ þ ϑðL β − L λ Þ gauge group, with divergences canceling between families. In particular, for the L μ − L τ model our result is in agreement with Refs. [59,75], Note that the finiteness of the one-loop result has important consequences for the leptophilic theories we consider. As pointed out in Ref. [51], the finiteness of ε implies that one is able to forbid tree-level kinetic mixing, albeit in a model dependent manner. This happens, for instance, when embedding the new leptophilic Uð1Þ group in a larger non-Abelian group G, which is completely independent from the SM sector. This choice of one-loop generated kinetic mixing should be seen as a conservative choice; in the absence of cancellation between tree and loop-level kinetic mixing, this yields the least constrained scenario for an L μ − L τ model. Additional constraints from first-family leptons are now relevant [72,76], especially ν − e scattering measurements, where the strength of the constraint makes up for the loop suppression in the coupling. For neutrino trident scattering, one can safely ignore loop-induced kinetic mixing contributions in the calculation since these are either smaller than the tree-level new physics contribution or yield very weak bounds compared to other processes. We emphasize that if accompanied by a consistent mechanism for the generation of the Z 0 mass terms and leptonic Yukawa terms, the models we consider constitute a UV complete extension of the SM. The treatment of such scenarios lies beyond the scope of this work, but we note that if their scalar sectors are light enough they can also yield rich phenomenology at low scales [77].

III. SIGNATURES OF LEPTONIC NEUTRAL CURRENTS
When a neutrino impinges on a detector, it has only two options for BSM scattering via a leptophilic mediator. In the simplest scenario, the neutrino interacts via the new mediator with the electrons of the detection medium. In this case, there is a tree-level ν − e scattering process that would be expected to show the clearest signs of new physics. For scattering off a hadron, however, the leptophilic nature of the mediator means that the first tree-level contribution will necessarily come from a diagram which also includes at least one additional SM mediator. Any neutrino-hadron scattering process can be embellished with the new boson to create a BSM signature. In general, the final states of these processes either will be identical to the original unembellished process (perhaps with missing energy) or will have an extra pair of leptons in the final state. These neutrino trilepton production processes, which we will refer to as tridents for simplicity, can be subdivided into four types: We note that these processes all occur in the SM, and so the hunt for new physics will necessarily be competing against a background of genuine SM events. Moreover, for final states with missing energy in the form of neutrinos, isolating a BSM signal would necessarily rely on spectral measurements, and other backgrounds have the potential to be large. In particular, the trident production of ννν and ννl will be seen as contributions to the neutral-current (NC) elastic and charged-current quasielastic (CCQE) processes, and we expect backgrounds to be insurmountable (see, e.g., Ref. [78] for new physics contributions to CCQE processes). The lll channels, on the other hand, are expected to have a much more manageable SM background. Trimuon production, for instance, has been measured in the past and provides a multitude of kinematical observables in the final state [79,80]. The SM rate for this channel contains radiative photon diagrams as well as hadronic contributions [81][82][83], while for leptophilic neutral bosons, the dominant contributions come from a weak process with initial and final state radiation of a Z 0 , making it a less sensitive probe of light new physics. Finally, the νll production, the most discussed trident signature in the literature, has already been observed in the dimuon channel [84][85][86]. This channel is by far the most important trident process for our study, as the leptonic subdiagrams contain only weak vertices in the SM.

A. Neutrino trident scattering
In the νll neutrino trident scattering, an initial neutrino scatters off a hadronic target producing a pair of charged leptons in the process. Since we focus solely on neutral current processes and on flavor conserving new physics, no mixed flavor tridents are relevant and we can write In the SM this process receives CC and NC contributions when α ¼ β and is a purely NC process if α ≠ β. The BSM contributions to trident production we consider are shown in Fig. 1. Beyond computing the Bethe-Heitler (BH) contributions considered previously, we show that radiative contributions to these processes are generally small. Using the narrow-width approximation (NWA), we compute the cross section for the radiation of a Z 0 particle from a neutrino-nucleus interaction, which can then promptly decay to an l þ l − pair. We call these contributions darkbremsstrahlung (DB) processes for their similarity with electron brehmsstrahlung in QED. We now discuss the two amplitudes individually. a. Bethe-Heitler. The BH amplitude can be written as follows: where Q 2 ≡ −q 2 ¼ ðP − P 0 Þ 2 is the momentum transfer and H μ EM the hadronic amplitude for coherent or diffractive electromagnetic scattering We refer the reader to Ref. [55] for the details on the treatment of the hadronic amplitude. The leptonic amplitude for NC scattering L μ reads FIG. 1. The BSM contributions to neutrino trident production considered in our calculation. The diagrams on the left are referred to as Bethe-Heitler contributions due to their resemblance to pair production. On the right, we show diagrams with a radiativelike Z 0 emission, which allows for the production of on-shell Z 0 particles, which subsequently decay into a charged-lepton pair.
In writing the equation above, we have introduced effective vector and axial couplings containing SM and BSM contributionŝ where A 's) are the SM vector (axial) couplings. Note the dependence on the positive kinematic variable K 2 in the BSM contribution, which can lead to a significant peaked behavior in the cross section. To avoid numerical difficulties, we have modified the phase space treatment proposed in [87,88], as shown in Appendix A.
b. Dark-bremsstrahlung. Because of the small decay width of the Z 0 (Γ ∝ g 02 M Z 0 ), one can obtain an estimate for its resonant production using the NWA. In the true narrowwidth limit, this process reduces to a 3-body phase space calculation and does not interfere with the BH amplitude. 4 where H μ W is the weak hadronic current (see Appendix B) and where ϵ Ã α ðk 2 Þ is the polarization vector of the Z 0 . The previous amplitude can then be squared and integrated over phase space for the total DB cross section. The different charged lepton final states can then be imposed with their respective branching ratios (BR). As a final remark, we note that the typical decay lengths of the new boson are typically below 1 cm for the parameter space of interest, such that their decay is indeed prompt.
From the previous discussions it is clear that the contributions to the total cross section at the lowest order in g 0 come from the interference between the BSM and the SM BH diagrams, and from the DB. The latter, however, contains an extra power of G F and is expected to be subdominant with respect to the BH interference. Our results for the individual flux integrated cross sections are shown in Fig. 2 for the μ þ μ − and in Fig. 3 for the e þ e − trident channels. We show the BH contributions as well as the DB one normalized by the SM trident cross section. All cross sections are flux integrated using the 62.4 GeV p þ DUNE flux described in Sec. IVA. For generality, we do not include the BR factors in the DB contributions, and so the green lines only apply for μ þ μ − tridents if M Z 0 > 2m μ and would suffer additional suppression due to the BR. In each figure we show two panels, one for vector couplings and one for axial-vector couplings. This is interesting from a purely computational point of view, as it shows explicitly the BH cross section scaling with the M Z 0 in the two cases. While the scaling is similar for dielectron tridents, it differs significantly between the vector and axial-vector cases of the dimuon cross section. This suggests the presence of mass suppression effects in the BH process. We do not investigate this further but note that there are large cancellations between the two BH diagrams in Fig. 1 which are only present for vectorlike couplings.
Finally, a cautionary remark on the axial-vector case. Despite the large enhancement present in the axial-vector case, we note that this is likely an artifact of our simplified model approach. In an UV completion, additional particles might contribute to the process, and these quadratic enhancements as a function of 1=M 2 Z 0 are expected to be regulated at some model dependent scale. It is beyond the scope of this paper to build such a model, and so for the sake of simplicity and concreteness, we only perform sensitivity studies for the vector model, where these enhancements are less problematic.

The equivalent photon approximation
We now comment on the EPA for neutrino trident production. This approximation is known to perform quite badly for the SM neutrino trident production cross section [55]. One may wonder, however, if the EPA gets better or worse when computing our BSM cross sections. Naturally, it would be most inadequate for the resonantlike cross sections, since the photon propagator and the strong 1=Q 4 behavior is absent. However, if one focuses on the BH contributions, a marginal improvement of the accuracy of the approximation is seen as one lowers the mass of the Z 0 mediator. In the SM, the ν − γ cross sections scale as a typical weak cross section, σ νγ ∝ G 2 Fŝ , whereŝ is the square of the center of mass energy of the ν − γ system. On the other hand, if the cross section is dominated by the BSM BH contributions, then as we take the limit of small Z 0 masses, it scales more similarly to a QED cross section, σ νγ ∝ 1=ŝ. This behavior, however, is only present at low masses and only for the BSM contribution. Since we are interested in regions of the parameter space where BSM and SM cross sections are of similar size, then we expect the total cross section to have a behavior that is a combination of the two. As a sanity check, we numerically verified that for parameter space points where the BSM contributions are of the same order as the SM cross section, the improvement in the accuracy of the EPA is still not satisfactory. For instance, the ratio between the EPA prediction and the full calculation for the dimuon channel assuming a Q max ¼ ð140 MeVÞ=A 1=3 goes from ≈30% in the SM to ≈60% for g 0 ¼ 8 × 10 −4 and M Z 0 ¼ 5 MeV. For this reason, we only use the full 2 → 4 calculation in what follows.

Trident kinematical distributions
The impact of new physics on the total cross section for trident production has been explored in the previous section. It is then natural to ask what the impact of new physics is on the kinematics of trident production which are, especially in the case of the invariant mass and angular variables, of utmost importance for background reduction. In this section we show how the new physics can alter the distributions of these important variables. All results that follow have been obtained using trident events produced by our dedicated Monte Carlo simulation (MC). Smearing and selection cuts have been applied as detailed in Sec. IV.
The variables of interest in background reduction are the charged lepton invariant mass m 2 ll ¼ ðp 3 þ p 4 Þ 2 and their separation angle Δθ ll . The invariant mass can be experimentally inferred from the energy of each charged-lepton and their separation angles, and so heavily relies on the experimental resolution to such parameters. In Fig. 4 we show the dimuon invariant mass spectrum between 4m 2 μ and 0.2 GeV 2 , and the dimuon separation angle between 2°a nd 18°for a light vector boson with M Z 0 ¼ 22 MeV. We show the results for the dielectron channel in Fig. 5. The light new physics here enhances these distributions at low values of these parameters. We show our results for two types of mediators, vector and axial-vector leptophilic bosons. Comparing the couplings necessary to produce similar BSM enhancements of the number of events, we see  that axial-vector bosons lead to larger enhancements with smaller couplings. In particular, it leads to greater spectral distortions for the Z 0 mass shown.

B. Neutrino-electron scattering
Neutrino-electron scattering has long been a valuable probe of both the SM and potential new physics [56,[90][91][92]. It is important to note that in the presence of novel leptophilic currents, experiments searching for e þ e − tridents would also observe anomalous ν − e event rates. In fact, given the larger statistics present in the ν − e scattering sample, this channel is expected to provide the leading constraints in our scenarios with tree-level couplings to electrons.
In order to compute the ν − e cross section in the presence of the new leptophilic interactions, we need to consider an analogous modification of the NC scattering amplitude where the vector (C V ) and axial (C A ) effective couplings include both the SM and the BSM contributions with, as usual, s W ≡ sin θ W , θ W being the weak angle, and T e is the kinetic energy of the recoil electron. The loop-induced kinetic mixing in the L μ − L τ model also induces a ν − e coupling The differential cross section is then given by where the left-and right-handed constants are given by For antineutrino scattering one obtains the cross section by exchanging C L α ↔ C R α . The kinetic energy of the outgoing electron is bounded by kinematics and the energy resolution of the detector, which effectively sets a threshold energy T th such that with T max ¼ 2E 2 ν =m e þ 2E ν , the maximum kinetic energy attainable. We define the effective total cross section for an initial neutrino energy E ν as This definition also ensures that the enhancement due to very light mediators becomes constant at around ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2m e T th p , as discussed in Ref. [56]. This is a consequence of the detector threshold and of the 2-body kinematics of the process. Finally, electroweak radiative corrections have been computed in the SM [93,94] but will not be included here. Since they correspond to a change of a few percent we do not expect them to affect our results very much.

ν − e kinematical distributions
The angle between the scattered electron and the neutrino beam direction, θ, is related to the electron energy as where y ≡ T e =E ν is the inelasticity (T th =E ν < y < 1) and E e ¼ T e þ m e is the outgoing electron energy. This implies that at OðGeVÞ neutrino energies, the electron recoil is very forward and obeys E e θ 2 < 2m e , up to detector resolution. For this reason, we choose to analyze our results in terms of E e θ 2 . In this case, the differential cross section becomes This distribution is particularly important for suppressing the background. Given the kinematics explained above, E e θ 2 must be smaller than 2m e for ν − e scattering, while it is often much larger for neutrino-nucleon scattering, the dominant background (see Sec. IVA). We show in Fig. 6 the expected ν − e event distribution as a function of E e θ 2 for the SM and a light Z 0 case, in the neutrino and antineutrino modes at the DUNE ND. As expected, the signal is extremely forward and the final distribution is highly sensitive to the angular resolution δθ of the detector. At a conservative value of δθ ¼ 1°, little information about the true distribution is left, and a significant portion of the signal lies in a region where E e θ 2 > 2m e . Therefore, shape information may improve the search for a light new physics only when the angular and energy resolutions of the detector are well understood.

C. Interference effects
Since for ν − e scattering and neutrino trident production there exists a SM contribution, we expect the experimental sensitivity to new physics to be dominated by the interference between SM and BSM contributions. We now argue what kind of interference one can expect in each one of these processes. For neutrino trident production we follow Ref. [95] and separate the differential cross section as where we dropped the flavor indices inV andÂ from (11b) for simplicity. This allows us to write the interference between the SM and the vector new physics as Depending on the region of phase space considered, the term proportional to dσ V−A can be of similar size to dσ V . However, dσ V−A changes sign as a function of the angular variables or energies, leading to small integrated cross sections (typically 2 orders of magnitude smaller than the integral of the dσ V term). Ignoring this term, one can then completely predict the type of interference in trident production. For ν μ → ν μ μ þ μ − trident production, for instance, C SM V > 0 and the second generation charge appears squared, leading to constructive interference in all cases. For ν μ → ν μ e þ e − trident events, on the other hand, C SM V < 0. If the first and second generation charges come in with opposite signs, then the interference is still constructive; otherwise destructive interference happens. The same considerations also apply to antineutrino scattering if one ignores the dσ V−A term. Finally, the axial-vector case is completely analogous taking V ↔ A in Eq. (22).
For ν − e scattering analytical expressions can easily be used [56]. Taking C SM Since y < 1, the interference term for ν μ − e is always positive (constructive), and forν μ − e it is always negative (destructive).

IV. DUNE SENSITIVITIES
Having studied the behavior of neutrino trident production and neutrino-electron scattering cross sections in the presence of light new bosons, we now apply our results in sensitivity studies for the DUNE ND. As discussed in Sec. II, we limit our studies to L e − L μ and L μ − L τ models with vector gauge bosons. We start with a discussion on the experimental details, highlighting the challenges of backgrounds and laying out our statistical methods in Sec. IVA. Then we show our main results in Secs. IV B and IV C, comparing our sensitivity curves to the leading bounds in the parameter space of the leptophilic models from other experiments.

A. Analysis techniques
The LBNF is expected to produce an intense beam of neutrinos and antineutrinos from a 1.2 MW proton beam colliding against a fixed target [29]. The DUNE ND, where the number of neutrino interactions is the largest, is expected to be located at a distance of 574 m from the target. Despite its design not being final yet [96,97], we focus on the possibility of a 75-t fiducial mass liquid argon (LAr) detector. Regarding the neutrino fluxes, we now concentrate on the option of a beam from 120 GeV protons with 1.1 × 10 21 protons on target (POT) per year. The LBNF could also provide higher or lower energy neutrinos depending on the proton energy, target, and focusing system used. We explore other possibilities shown in Table I, and we take the flux files provided in Refs. [98,99]. We assume that the experiment will run 5 years in neutrino mode and another 5 years in antineutrino mode. The final exposure, therefore, will vary with beam designs and is equal to a total of 11 × 10 22 POT in the case of 120 GeV protons. To generate neutrino scattering events, we use our own dedicated MC simulation, Gaussian smearing the true MC simulation energies and angles as a proxy for the detector effects during reconstruction. We assume an energy resolution of ffiffiffi ffi E p Þ for e=γ showers (muons) and angular resolutions of δθ ¼ 1°for all particles [100].
An interesting addition to the design of the DUNE ND would be a magnetized high-pressure GAr tracker placed directly behind the LAr module [101]. The lower thresholds for particle reconstruction and the presence of a magnetic field is expected to improve event reconstruction and reduce backgrounds to neutrino-electron scattering and neutrino trident production. We note that despite the relatively small fiducial mass of such a GAr module, ≲1 tonne, it would still provide a sizable number of these rare leptonic neutrino scattering processes. With the intense flux at DUNE and the large number of POT, the ν − e scattering measurement will not be statistically limited, with order 10 4 events in the DUNE ND after a few years. Systematics from the beam and detector are then the limiting factor for the sensitivity to new physics in this measurement. Current work on neutrino flux uncertainties shows that normalization uncertainties can be reduced to the order of 5% [102][103][104], with similar projections for DUNE [29]. The electron energy threshold also plays a role in the new physics search. In particular, for new light bosons the enhancement at very low momentum transfer 2T e m e has a cutoff at the minimum electron recoil energy [see Eq. (19)]. This implies that the experiment is no longer sensitive to the Z 0 mass below ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2T th m e p . In our analysis, we assume a realistic overall normalization systematic uncertainty of 5% and a ν − e scattering electron kinetic energy threshold of 600 MeV.
Lowering systematic uncertainties on the flux is challenging given the large hadroproduction and focusing uncertainties at the LBNF beam. Here, improvements on the experimental side in determining the neutrino flux will be extremely valuable (see, e.g., Ref. [105]). If one is searching for novel leptophilic neutral currents, hadronic processes and inverse muon decay measurements are available, but these are limited either by theoretical uncertainties or by statistics and might not be applicable in the whole energy region of interest. As to the electron energy, assuming a threshold as low as 30 MeV would be safe for electron detection, but at these low energies backgrounds can be incredibly challenging due to the overwhelming π 0 backgrounds. Increasing this threshold to 600 MeV, however, has little impact on our sensitivities and is only 200 MeV below the threshold used in the most recent MINERνA analysis [53], where good reconstruction is important for measuring the flux. For e þ e − and μ þ μ − tridents, we refrain from increasing the analysis thresholds from a naive 30 MeV. This is certainly an aggressive assumption but it is necessary if e þ e − tridents are to be measured, since these events are quite soft [55]. Thresholds for μ þ μ − tridents are much less important since the events are generally more energetic than their dielectron analogue.
a. Backgrounds (ν μ → ν μ l þ l − ).-We now discuss the individual sources of backgrounds to neutrino trident production. A pair of charged leptons is very rarely produced in neutrino interactions, usually coming from heavy resonance decays [24,86,[106][107][108]. Since our signal is mostly coming from coherent interactions with nuclei, cuts in the hadronic energy deposition in the detector E had , often large in heavy meson production processes, can help reduce backgrounds. Coherent and diffractive production of mesons is an exception to this, in particular pion production [109][110][111][112], which is the main background to trident due to particle misidentification (misID). Muons are known to be easily spoofed by charged pions, making CC ν μ interactions with π AE in the final state (CC1π) one of the largest contributions to the backgrounds of μ þ μ − tridents. Similarly, NC π 0 production stands as the leading background to e þ e − tridents when the photons are misIDed as two electrons, or if one of the photons pair converts and the other escapes detection. In Ref. [55], we have shown that the μ þ μ − and e þ e − pairs produced in trident have small separation angles (Δθ) and possess small invariant masses (m 2 ll ), and that both charged leptons are produced with small angles with respect to the neutrino beam (θ AE ). With simplified misID rates, we used the GENIE [113] event generator to show that simple kinematical cuts can reduce backgrounds significantly, achieving a significance of S μμ = ffiffiffiffiffiffiffi B μμ p ∼ 44 and S ee = ffiffiffiffiffiffiffi B ee p ∼ 17.3 for the DUNE ND in neutrino mode, where S and B stand for signal and background, respectively. In our current analysis we implement the same kinematical cuts, which are as follows: m 2 μμ < 0.2 GeV 2 , θ AE < 15°, and Δθ < 20°for the μ þ μ − channel; and m 2 ee < 0.1 GeV 2 , θ AE < 20°, and Δθ < 40°for the e þ e − one. We impose these cuts again in our signal analysis and point out that the new physics enhancement happens precisely in this favorable kinematical region (see Sec. III A 2). The degree with which the experiment will be able to reduce backgrounds will rely on reconstruction properties of the signal and background final states. In particular, the detector containment of the charged-lepton pairs, as well as pions and photons, is crucial for momentum and invariant mass reconstruction, and so a detector simulation is desirable. Since we do not aim to develop a full experimental analysis and since the DUNE ND design is still under debate, we present our results with no backgrounds in Fig. 9 and vary the total background rate in Fig. 10, all the while applying the cuts above. This illustrates the impact of worse detector performance in background rejection.
b. Backgrounds (ν − e).-For neutrino-electron scattering, backgrounds will arise from either the genuine production of an electron or via the misID of particle showers in the detector, both in the absence of observable hadronic energy deposition. The former scenario happens mostly by the CC interactions of the flux suppressed ν e states present in the beam. The main contribution will be from CCQE interactions where the struck nucleon is invisible either for being below threshold or due to nuclear reabsorption. The misID of a photon initiated EM shower for an electron one is expected to be rare in LAr, where the first few centimeters of the showers can be used to separate electrons and photons by their characteristic dE=dx. However, the large NC rates for the production of single photons and π 0 can become a non-negligible background.
For instance, coherent NC π 0 production leaves no observable hadronic signature and may look like a single electron if one of the photons is misidentified and the other escapes detection. Finally, after misID happens, the signal can still look unique in its kinematical properties. In particular, E e θ 2 cuts can dramatically reduce backgrounds due to the forwardness of our signal (see, e.g., [53,114]).
c. Statistics.-In order to assess the potential of DUNE to discover new physics, we perform a sensitivity analysis using a χ 2 test with a pull method for systematic uncertainties. Our goal is to assess when DUNE would be able to rule out the SM, and so we generate BSM events and fit the SM prediction to it. Our χ 2 function is defined as where the number of events for the BSM case is given by N BSM , the SM number of events is N SM , and the number of background events is N BKG . The nuisance parameters α and β, with their uncertainties σ norm and σ BKG , take into account normalization uncertainties from the flux and detector, and uncertainties on the background prediction, respectively. For the DUNE ND, we assume σ norm ¼ 5% and σ BKG ¼ 10%. These systematics will likely be dominated by flux normalization uncertainties and can only be measured with interactions that do not depend on the leptophilic BSM physics.
New vector bosons with couplings to the first and second generation leptons can be probed very effectively in neutrino experiments by measuring the ν − e scattering rate. This has been recognized in the literature [51,56,115], where bounds from various experiments, including CHARM-II [116], TEXONO [117][118][119], and Borexino [120] have been derived on these bosons. Curiously, the bound calculated from the CHARM-II data has been pointed out by Ref. [51] to be too optimistic. The uncertainty on the neutrino flux is a real hindrance for these measurements which has not been taken into account when these bounds were computed. This is particularly important for measurements with large statistics, and for this reason we do not show the CHARM-II bound here. The measurement ofν e − e scattering at TEXONO, on the other hand, is statistically limited, and the bound it places on this class of models can safely ignore the flux systematics. This turns out to provide the strongest limit in a large region of the L e − L μ parameter space. Trident bounds can be obtained for this model, but due to their lower statistics and more involved kinematics, are subdominant.
We show our results for the DUNE ND in Fig. 7. Our results are for the combined ν þν modes and do not include backgrounds. The opposite charges between the first and second families imply constructive interference between the SM and BSM contributions for neutrino scattering, contrary to what happens in a B − L model, for instance. Therefore, the strongest bounds on this model can be obtained at DUNE in neutrino mode. It is clear, however, that the degree with which DUNE can probe unexplored parameter space is a question of how much the  [117][118][119], searches at the BABAR e þ e − collider [121,122], and beam dump experiments [51]. uncertainties on the flux can be lowered. To illustrate this effect, we vary the normalization systematics on the right panel of Fig. 8, going from a conservative 10% to an aggressive 1% uncertainty. The effect of changing the thresholds is very small, being most important in a region already probed by other experiments. Different beam designs seem to have only a small impact on the sensitivity, as shown on the left panel of Fig. 8.
Since we show the bounds obtained from the neutrino and antineutrino runs combined, it is not possible to see the effects of destructive interference. If only channels with destructive interference were available, however, it would have been possible to allow for cancellations between the total interference and the square of the BSM contributions in certain regions of parameter space at the level of the total rate. The region where this cancellation happens depends strongly on the neutrino energies involved and on the integrated phase space of the recoiled electron. In that case, one expects that the sensitivity to the lowest new physics couplings comes, in fact, from the search for a deficit of ν − e scattering events, as opposed to the constructive interference case where an excess of events is always produced. We note that this has no significant impact on the sensitivity of a leptophilic Z 0 , but might provide crucial information about the nature of the Z 0 charges in case of detection.
The trident bounds we obtain are not competitive for this model despite the fact that the trident cross sections receive similar enhancements to that of ν − e scattering. This is due to two reasons: the low number of events and the nontrivial kinematics of trident processes. Since the neutrino is essentially scattering off virtual charged leptons produced in the Coulomb field of the nucleus, it has to typically transfer more energy to the system than it would in a scattering off real particles in order to produce visible signatures. This remark also helps us to explain the behavior of the sensitivity curves at the lowest masses. While ν − e scattering cross sections become insensitive to the boson mass at ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2m e T th p , the trident cross sections do not. This behavior is most dramatic in the e þ e − tridents, but is also present in the μ þ μ − one. This is a consequence of the 4-body phase space kinematics, where now the momentum transfer through the Z 0 propagator is no longer trivially related to the final state particle energies, as in 2 → 2 processes. It should be noted, however, that both the dimuon and the dielectron trident rates become nearly independent of M Z 0 below the muon and the electron mass, respectively, where only a logarithmic dependence is expected [52].
DUNE can also probe this class of models in a different way. In the context of long range forces in neutrino oscillation experiments and with the same choice of charges, Ref. [123] places competitive bounds in this model with Super-Kamiokande data and makes projections for DUNE. The matter potential created by the local matter density modifies the dispersion relation of the neutrinos with lepton nonuniversal charges, leading to very competitive bounds in our region of interest. Similar considerations have also been explored in the context of high-energy astrophysical neutrinos [124]. Other experimental searches have been conducted at electron beam dumps. This technique consists of producing the Z 0 boson at the target via radiative processes such as e þ A → e þ A þ Z 0 , and we look for the visible decays of the boson in the detector. In this model, the decay products are mostly e þ e − states and the bounds are only applicable at appreciably small values of g 0 and M Z 0 , where the lifetime of the Z 0 is sufficiently large. Probing the large mass region, on the other hand, requires high-energy experiments. In that regime, the strongest bounds come from searches at the e þ e − collider BABAR. These come about in two FIG. 8. The ν − e scattering sensitivity to the L e − L μ model at 90% C.L. In the left panel we show the sensitivity using different choices for the neutrino flux, and in the right we use the neutrino beam from 120 GeV protons and vary the normalization systematic uncertainty from an aggressive 1% to a conservative 10%. ways: looking for the visible decay products of a Z 0 produced radiatively or in heavy meson decays [121], or exploring the BR into invisible final states [122].
In this section we evaluate the DUNE ND sensitivity to the presence of a light vector Z 0 charged under L μ − L τ . Beyond being anomaly free, this choice of charges allows for positive contributions to the anomalous magnetic moment of the muon, a μ ¼ ðg − 2Þ μ , as discussed in Refs. [72,76,[125][126][127]. This quantity is well known for a ∼3.7σ discrepancy between the experimental measurement [128] and the theory predictions [129,130]. If future efforts to measure it [131] confirm this disagreement and if theoretical uncertainties are better controlled in the next few years, then constraining new physics scenarios that could contribute to a μ is of utmost importance.
This model can significantly impact neutrino trident production of a muon pair. In fact, the leading bound in this parameter space for masses M Z 0 ≲ 200 MeV comes from the CCFR measurement of the same neutrino trident channel [85]. CCFR observed 37.0 AE 12.4 events, extracting a measurement of the trident cross section of σ CCFR =σ SM ¼ 0.82 AE 0.28. Curiously, the measurement by CHARM-II [84] provides weaker constraints on this model despite seeing a larger number of trident events, namely 55 AE 16 events in total, most likely due to the 1σ upward fluctuation of the measurement: σ CHARM−II =σ SM ¼ 1.58 AE 0.57. Other important bounds from ν − e scattering have also been obtained using the kinetic mixing parameter generated at one loop, the strongest of which uses data from Borexino [58] and is relevant only for the low mass region M Z 0 ≲ 20 MeV.
At DUNE, both of these measurements are possible, allowing one to constrain this model in different ways.
FIG. 9. The DUNE ND neutrino scattering sensitivities for L μ − L τ at 90% C.L. The upper panel shows the case with no kinetic mixing, and the lower panel shows the case with the loop-induced mixing. Bounds from neutrino-electron scattering apply only to the latter. We also show bounds from BABAR [132], from LHC [133], from Borexino [58], and from the neutrino trident production measurement at CCFR [52,85]. Recent cosmological bounds for the two kinetic mixing cases derived in Ref. [75] are also shown.
We show our results in Fig. 9, without including backgrounds. In this scenario, DUNE would be able to cover all the 2σ region compatible with the ðg − 2Þ μ measurement only with the μ þ μ − trident events. For the low mass region, measuring the ν − e scattering rate can provide a complementary probe of this region, depending most strongly on the systematic uncertainties DUNE can achieve. We note that analysis thresholds used for ν − e scattering have little impact on the sensitivity in the region of interest. Our conclusion that DUNE can cover all of the ðg − 2Þ μ region holds provided backgrounds are kept below the SM signal rate. This can be seen when we include backgrounds with different assumptions in the right panel of Fig. 10. Finally, different assumptions for the beam design have little impact on the sensitivity, as show in the left panel of Fig. 10.
Apart from neutrino scattering, dedicated searches for resonances decaying into μ þ μ − in four muon final states have been performed at BABAR [132], looking for e þ e − → μ þ μ − Z 0 ð→ μ þ μ − Þ. At the LHC, the Z → 4μ measurement performed by the ATLAS Collaboration [133] was used to derive a constraint in the L μ − L τ parameter space in Ref. [52]. Recently, the CMS Collaboration performed a dedicated search for a resonance between M Z 0 ¼ 5 and 70 GeV, significantly improving previous constraints at large masses [134]. Big bang nucleosynthesis bounds were studied in [72,127] and shown to constrain the mass of the boson to be M Z 0 ≳ 5 MeV. Recently, additional constraints from cosmology were derived given that the presence of very light Z 0 bosons changes the evolution of the early universe [75]. In particular, the decays and inverse decays induced by the new leptophilic interactions can modify the neutrino relativistic degrees of freedom, requiring M Z 0 ≳ 10 MeV in order for ΔN eff < 0.5 for the case with no kinetic mixing. The authors of Ref. [75] also found that an additional Z 0 boson can alleviate the tension in the different measurements of the Hubble parameter. Let us stress here that all these bounds will be complementary to possible future constraints that can be obtained by the DUNE program, as shown in Fig. 9.

V. CONCLUSIONS
Although the next generation neutrino oscillation experiments are primarily designed for making precision measurements of the neutrino mixing parameters, the unprecedented fluxes and large detectors will allow for many nonminimal new physics searches. In this work, we have considered the physics potential of the DUNE ND for constraining the existence of an additional anomaly-free Uð1Þ gauge group giving rise to a Z 0 boson which only couples to leptons-a form of a purely leptophilic neutral current. Specifically, we have considered the anomaly-free scenarios with charges associated with the lepton number difference L α − L β . Focusing on the two most promising neutrino scattering processes, ν − e and νll trident scattering, we have computed expected sensitivity curves for the DUNE ND for a variety of charge assignments.
In performing our sensitivity studies as a function of the coupling and mass of the Z 0 boson, we have remained as faithful as possible to the real experimental conditions of a LAr detector. Our main results rely on the realistic assumptions of flux uncertainties of 5% and feasible exposures. To avoid large backgrounds, we have also implemented kinematical cuts on the neutrino trident sample, and a kinetic energy threshold of 600 MeV for ν − e scattering events. The parameter space that can be probed by ν − e scattering in the L e − L μ scenario is at least 2 times better than the e þ e − and almost 20 times better than the μ þ μ − trident channels, especially for the lower mass region. In this case, the DUNE ND would improve only FIG. 10. The dimuon neutrino trident sensitivity to the L μ − L τ model with no kinetic mixing at 90% C.L. In the left panel we show the sensitivity using different choices for the neutrino flux, and in the right we use the neutrino beam from 120 GeV protons and scale the background with respect to the total number of SM trident events after cuts. slightly on previous ν − e scattering bounds, especially at around M Z 0 ∼ 100 MeV. We do not expect e þ e − trident measurements at DUNE to improve our coverage of the L e − L μ Z 0 parameter space, but note this process has a distinct dependence on M Z 0 if compared to ν − e scattering.
If the light vector Z 0 is charged under L μ − L τ , we have found that the dimuon trident measurement could provide the leading bound in this parameter space. This is particularly interesting as these models can also explain the discrepancy between the measurement of the anomalous magnetic moment of the muon and its SM prediction. We expect that DUNE will be able to fully explore the ðg − 2Þ μ motivated parameter space provided backgrounds are kept under control. The robustness of our results is tested against different choices of neutrino fluxes, where we find that despite the larger rates at higher neutrino energies and the larger BSM enhancement at lower energies, the sensitivities are very similar.
Improvements to the experimental sensitivities we have displayed in Figs. 7 and 9 can be achieved by reducing uncertainties on the neutrino flux and detection. From the experimental side, novel detection techniques suitable to rare neutrino events are currently under discussion, such as the magnetized HPgTPC [101] and the Straw Tube Tracker [135,136]. Together with improved analysis techniques, these will help to improve upon our projections for the sensitivity of DUNE to new physics that might be hiding at light masses and small couplings.
where Q W ¼ ð1 − 4s 2 w ÞZ − N and F W ðQ 2 Þ stands for the weak form factor of the nucleus. We implement the Helm form factor as in [138], defined as where j 1 ðxÞ stands for the spherical Bessel function of the first kind, s ¼ 0.9 fm and R ¼ 3.9 fm for 40 Ar.