Search for a light muon-philic $Z^\prime$ with the NA64-$e$ experiment at CERN

The extension of Standard Model made by inclusion of additional $U(1)$ gauge $L_\mu-L_\tau$ symmetry can explain the difference between the measured and the predicted value of the muon magnetic moment and solve the tension in $B$ meson decays. This model predicts the existence of a new, light $Z^\prime$ vector boson, predominantly coupled to second and third generation leptons, whose interaction with electrons is due to a loop mechanism involving muons and taus. In this work, we present a rigorous evaluation of the upper limits in the $Z^\prime$ parameter space, obtained from the analysis of the data collected by the NA64-$e$ experiment at CERN SPS, that performed a search for light dark matter with $2.84\times10^{11}$ electrons impinging with 100 GeV on an active thick target. The resulting limits, despite being included in a region already investigated by neutrino experiments,touch the muon $g-2$ preferred band for values of the $Z^\prime$ mass of order of 1 MeV. The sensitivity projections for the future high-statistics NA64-$e$ runs demonstrate the power of the electrons/positron beam approach in this theoretical scenario.


INTRODUCTION
The Standard Model (SM) of particle physics works remarkably well in describing and interpreting the experimental results provided by different, complementary efforts, operating at different energy scales [1]. However, oretical prediction computed by the Muon g − 2 Theory Initiative [4], a µ (Exp) − a µ (SM) = (251 ± 59) × 10 −11 . Even if recent alternative results obtained through lattice QCD calculations may possibly release this tension [5], still the muon magnetic moment anomaly motivates the interest toward beyond SM scenarios that may explain it. In this work, we consider the SM extension in which the anomaly-free combination L µ − L τ is associated to a new U (1) gauge symmetry, thus introducing a new massive vector boson Z coupled to the difference between the second and third generation leptonic currents [6,7]. The corresponding new lagrangian terms read [8]: where Z µν ≡ ∂ µ Z ν − ∂ ν Z µ is the Z field strength, m Z is the Z mass, P L = (1 − γ 5 )/2, and g Z is the coupling between the Z boson and the L µ − L τ SM current. At leading order, the Z contributes to the muon magnetic moment as: where F (x) = 1 0 dz 2z(1−z) 2 (1−z) 2 +x 2 z [9][10][11]. This can explain the observed muon magnetic moment discrepancy if the parameter m Z and g Z lie in a well defined area of the parameters space, roughly defined by g Z ∈ [3.2, 5.5] × 10 −4 at 2σ for m Z m µ [12]. The lack of a tree-level coupling with the SM electron also means that the model does not contribute appreciably to the e − magnetic moment, in agreement with the experimental observation [13]. The Z model, either in the "vanilla" form described before or in association with more elaborated SM extensions, has also been advocated to explain other SM anomalies, such as the B decay anomaly [14][15][16][17] and the lepton-flavor universality violation [18][19][20]. The Z model has also been connected to the Dark Matter (DM) [10,[21][22][23][24][25][26] and to the neutrino mass phenomenology [27].
These arguments recently motivated a large number of complementary efforts to search for the Z , either by performing a re-analysis of existing experimental datasets, or proposing new, dedicated experiments. The BaBar [28] and CMS [29] experiments investigated existence of the Z by exploiting the visible decay channel Z → µ + µ − , searching for a resonance peak in the dimuon mass distribution, on top of the SM background. The Belle-II experiment focused instead on the Z invisible decay channel, exploiting the reaction e + e − → µ + µ − Z , where the Z is radiated by one of the final state muons, searching for a resonance peak in the recoil mass of the final state muons [30]. Dedicated Z searches at Belle-II exploiting mono-photon signatures have also been suggested [31][32][33]. Stringent upper limits on the g Z coupling have also been obtained by neutrino experiments, such as CCFR [34], Borexino [35,36], and COHERENT [37]. Among future proposals, FASER-ν at CERN aims to search for Z via neutral-current deep-inelastic neutrinonucleon scattering [38].
Since the Z couples predominantly to second and third generation leptons, the most effective experimental strategy to investigate this model at accelerators is by exploiting muon beams. Dedicated efforts have been proposed at CERN (NA64-µ [39][40][41][42][43]) and Fermilab (FNAL-µ [40] and M 3 [25]), with NA64-µ having already completed a pilot run in October 2021. Nevertheless, thanks to the presence of a loop-induced Z -electron coupling, e − beam experiments can also probe a significant portion of the Z parameter space, somehow paving the road to next generation efforts.
In this work, we present the upper limits introduced to the Z parameters space from a search performed with the NA64-e experiment at CERN [44], in the mass range 1 − 600 MeV. We consider two Z models, the "vanilla" one and a "dark" scenario in which the Z couples predominantly to light dark sector particles. In the second case, similarly to what was done in Refs. [22,25,45], we introduce a dark scalar particle with mass m χ and coupling to Z defined by the following Lagrangian: where J µ D is the dark vector current given by: We assume the mass hierarchy m χ < m Z /2 and the couplings ratio g D /g Z 1. This choice results to a preferred combination of the model parameters that can reproduce the DM relic density observed at present, in the hypothesis that χ particles are responsible for it: where f = 1-10 depends on the nature of the dark sector particle (scalar, Dirac or Majorana fermion, . . .) [25]. We observe that, since the L µ − L τ model considered here is associated to a U (1) gauge symmetry, all the Z interactions should be proportional to the gauge coupling g Z . Hence, g D should be decomposed into a product g Z · q χ , where q χ is the χ field charge associated to the new U (1) gauge group. The "dark" scenario considered in this work corresponds to the case q χ 1, already introduced in Ref. [25]. For convenience, all the results reported below will be presented as a function of g D .
The paper is organized as follows. In Sec. II we introduce the phenomenology for Z production in fixed target electron/positron beam experiments, with a focus on missing-energy efforts. In Sec. III we briefly describe the NA64-e experiment at CERN. In Sec. IV we present our strategy to extend the existing NA64-e results to the Z model, for both the "vanilla" case and the "dark" one, and finally in Sec. V we show and discuss the obtained results, including the sensitivity projections for future NA64-e runs.

II. Z PRODUCTION IN FIXED TARGET ELECTRON-BEAM EXPERIMENTS
The Z model considered in this work does not include explicitly a kinetic mixing term between the Z and the SM photon, that would result in a tree-level coupling with the SM electric charge. However, such a coupling arises naturally from the one-loop diagrams reported in Fig. 1, introducing an effective e ± − Z interaction term eΠ(q 2 )Z µ (eγ µ e), where the complex function Π(q 2 ) depends on the momentum q 2 carried by the Z [45][46][47]: A plot of the Π(q 2 ) function is reported in Fig. 2. As already pointed out in Ref. [46], the dependence of the e ± − Z coupling on the momentum is a unique feature of this model, and makes the phenomenology of Z searches at electron/positron beam experiments significantly different than the dark photon case, where the coupling is constant [48]. For very small values of the Z momentum, q 2 m 2 µ , the Π(q 2 ) function assumes the constant value Π(0) = e g Z 6π 2 ln mτ mµ 0.0144 · g Z , and the Z − e ± interaction resembles that of the "traditional" dark photon model under the exchange ε ↔ 0.0144 · g Z , where ε is the dark photon kinetic mixing parameter [48,49]. At larger momentum values, however, there is an enhancement of Π(q 2 ), with a maximum value for q 2 = 4m 2 µ , where its magnitude is a factor ∼ 1.5 larger than its small-momentum value, resulting in a significant increase of the Z production yield in this kinematic region.
The main Z production processes in the collision of a high energy electron or positron beam with a fixed thin target are shown in Fig. 3. Diagram (a) corresponds to the so-called Z -strahlung process, in which a Z is radiatively emitted by the lepton interacting with the electromagnetic field of a nucleus in the target. Diagrams (b) and (c), relevant only for an impinging positron, correspond to the non-resonant (b) and resonant (c) e + e − annihilation. These processes are analogous to those relevant for the production of a dark photon, with the significant difference of the loop-induced Π(q 2 ) factor appearing in the effective coupling. The corresponding production cross-section σ P (E, m Z ) formulas, con- sidering an on-shell Z and a beam energy E, can be obtained from the corresponding expressions for a dark photon model (see, e.g., [50][51][52]) with the substitution ε ↔ Π(m 2 Z ). Depending on the model and on the specific parameter values, the produced Z can decay to different final states. For the "dark" case, in the mass range 2m µ < m Z < 2m τ , the following decay channels are possible (neglecting the strongly suppressed Z → e + e − decay): where α Z ≡ , r µ ≡ m µ /m Z , r χ ≡ m χ /m Z , and the neutrino channel refers to the summed contributions from ν µ and ν τ . The "vanilla" scenario results can be simply obtained by setting g D = 0. We observe that, for m Z < 2m µ , only the invisible decay channels to neutrinos, and eventually to dark sector particles, are allowed.
For a thin target electron-beam experiment, with t X 0 , where t is the target thickness and X 0 the radiation length, only diagram (a) contributes to the Z event yield, scaling as: where x denotes a set of kinematic variables to describe the final state phase space.
In the case of a thick target electron-beam experiment with t X 0 , all the aforementioned production channels contribute to the Z yield due to the presence of the to the χ − e scatterings in the detector) was computed as: where N χχ is the total number of LDM particles (χ + χ) propagating from the beam-dump and impinging on the detector, L det and n e (N Av /A ρZ) are the detector length and the electron density, respectively, ε s is the average signal detection efficiency, and σ * χe is the total χ − e scattering cross-section integrated over recoil electron energies larger than the detection threshold E thr . N χχ was computed by projecting the χ angular distribution in the dump to the detector front-face plane and measuring the fraction of crossing particles. To evaluate σ * χe and to determine the energy and angular spectrum of recoiling electrons, we used the differential cross-section reported in Ref. [7]: The detection efficiency ε s of each experim sidered was determined by applying the sa cuts used in the original analyses. Furthe given in the following.
E137 is a BDE that ran at SLAC in 1980ing for long-lived neutral objects which m duced in the electromagnetic shower initiate electrons in the SLAC Beam Dump East. T rameters of the experiment are summarize secondary electrons and positrons in the electromagnetic shower induced by the primary electron. For the specific case of a missing-energy experiment, in which the signal is associated to the production of invisibly-decaying Z particles with energy greater than a threshold E M iss cut , the event yield scales as [51,53]: where T − (T + ) is the secondary electrons (positrons) differential track-length distribution [54,55] as a function of their energy E, and T ± ≡ T − + T + . The quantity is the differential cross section per nucleus for radiative Z production with respect to the final state invisible particles total energy E F , while σ Res is the total cross section per electron for the resonant production process. Z is the atomic number of the target material 1 . The two integrals in the radiative contribution term are performed over the E F > E M iss cut range. For the resonant production the kinematic constraint E F = E reduces the dimensions of the integral region. We did not include the non-resonant Z production mechanism in this computation, since for a primary electron beam the contribution to the total yield due to the annihilation with secondary positrons (see e.g. Fig. 3(b)) is negligible [51].

III. THE NA64-e EXPERIMENT
The NA64-e experiment at CERN is devoted to the search for dark sector particles feebly interacting with electrons. NA64-e exploits the 100 GeV high-purity, lowcurrent electron beam from the H4 beamline to perform the search, by measuring event-by-event the energy deposited in a thick active target, looking for events with large missing energy (see Refs. [44,[56][57][58] for a complete description of the detector and of the missing-energy approach). A schematic view of the detector is shown in Fig. 4. The NA64-e active thick target is an inhomogeneous electromagnetic calorimeter (ECAL), with energy resolution σ E /E 10%/ E(GeV) + 4%. The ECAL is made by 150 alternated layers of 1.5-mm thick lead plates and 1.5-mm thick plastic scintillator tiles, for a total length of about 40 X 0 ; the detector is segmented into a 6x6 matrix of independent transverse cells, with each cell further divided into a 4 X 0 pre-shower section and a main section.
Two main types of backgrounds, resulting to missing energy events, affect the NA64-e experiment. The first type is associated to the production of one or more penetrating particles in the ECAL by the primary beam. To suppress this contribution, a massive hadronic calorimeter (HCAL) is installed downstream of the ECAL. A high-efficiency plastic scintillator detector (VETO) is also located between the two calorimeters, to identify events in which penetrating charged particles are produced. The second type of background events is due to residual ≈ 1% hadron contaminants in the primary beam. To suppress these, a syncrotron radiation beam-tagging system (SRD) is installed upstream of the ECAL [59]. The NA64-e detector assembly also includes a magnetic spectrometer to measure the momentum of impinging particles. This consists of two successive dipole magnets (total magnetic strength Bdl 7 T· m) and a set of upstream and downstream tracking detectors, Micromegas (MM), Strawtubes (ST) and Gaseous Electron Multipliers (GEM). Finally, a set of beam-defining plastic-scintillator counters (SC) is present. During operations, the majority of the primary electrons gives rise to Compared to the analysis of Ref. [38], a number of improvements, in particular, in the track reconstruction were made in the 2018 run to increase the overall efficiency. Also, the zero-degree calorimeter HCAL 0 was used to reject events accompanied by hard neutrals from the upstream e − interactions; see Fig. 1.
In order to avoid biases in the determination of selection criteria for signal events, a blind analysis was performed. Candidate events were requested to have the missing GeV; E HCAL < 1 GeV) was defined based on the energy spectrum calculations for A 0 s emitted by e AE from the electromagnetic (e-m) shower generated by the primary e − s in the target [48,49]. A Geant4 [50,51] based Monte Carlo (MC) simulation used to study the detector performance, signal acceptance, and background level, as well as the analysis procedure including selection of cuts and estimate of the sensitivity are described in detail in Ref. [38].
The left panel in Fig. 2 shows the distribution of ≃3 × 10 4 events from the reaction e − Z → anything in the ðE ECAL ; E HCAL Þ plane measured with loose selection criteria requiring mainly the presence of a beam e − identified with the SR tag. Events from area I originate from the QED dimuon production, dominated by the reaction e − Z → e − Zγ; γ → μ þ μ − with a hard bremsstrahlung photon conversion on a target nucleus and characterized by the energy of ≃10 GeV deposited by the dimuon pair in the HCAL. This rare process was used as a benchmark allowing us to verify the reliability of the MC simulation, correct the signal acceptance, cross-check systematic uncertainties, and background estimate [38]. Region II shows the SM events from the hadron electroproduction in the target that satisfy the energy conservation E ECAL þ E HCAL ≃ 100 GeV within the energy resolution of the detectors.
Finally, the following selection criteria were chosen to maximize the acceptance for signal events and to minimize background. an electromagnetic shower in the ECAL, with full energy release in the detector. To suppress the rate of events processed and written to disk, the experiment trigger requires, other than a coincidence signal between the SCs, that the total energy sum signal from the ECAL corresponds to an energy deposition of less than 80 GeV.
NA64-e already completed data-taking campaigns in 2016, 2017, 2018 with a total accumulated charge of about N EOT = 2.84 · 10 11 electrons-on-target (EOT). The selection criteria adopted in the analysis were identified by optimizing the experiment sensitivity, adopting a blind-analysis approach [44,58,60]. These include the requirement to have a well reconstructed track in the upstream spectrometer, with momentum in the range (100 ± 3) GeV, an in-time cluster in the SRD detector, and a shower signal in the ECAL with the longitudinal and transverse shape of a missing-energy event. The latter selection also included a 0.5 GeV energy cut for the ECAL pre-shower section. After applying all selection cuts, no events were observed in the signal region, defined by the two requirements E ECAL < 50 GeV and E HCAL < 1 GeV; this observation is compatible with the estimate of (0.53 ± 0.17) background events, mostly due to the interaction of electrons with upstream beamline elements, producing a soft electron hitting the ECAL and one or more hadrons at large angle missing the NA64-e detector. This result was used to set an exclusion limit for the production of an invisibly decaying dark photon (A ), taking into account both the radiative and resonant production [60].

IV. METHODOLOGY
The analysis presented in this work is based on the dataset already scrutinized by the NA64 collaboration to set limits for an invisible decaying dark photon, preventing us to adopt a blind approach. For this work, we decided to follow a strategy based on the same selection criteria adopted in the aforementioned analysis, including the signal region definition. Our approach is based on the observation that, for an electron beam missing energy experiment, the Z signal would differ from the A one only due to the momentum dependence of the e − e − Z vertex, associated to the Π(q 2 ) function.
We considered first the simpler case in which only the radiative Z production channel is included. The 90% C.L. upper limit ε up for the dark photon kinetic mixing parameter reported by NA64-e in Ref. [44], given the negligible number of expected background events, corresponds to an expected number of signal events equal to N up 2.3 via the relation: where dσ A ,Rad dE F is the differential A production cross section per nucleus divided by ε 2 , η A is the corresponding signal acceptance and detection efficiency, and N is the overall normalization factor, accounting for the total accumulated number of EOT and for the detector material composition. Due to the detector geometry, T (E) is almost the same for the plastic scintillator and the lead, while a Z 2 dependence is included in the cross section: since N scales as the material density over the atomic mass, in the following we will consider only the contribution from the lead (the same approximation was made in the analysis presented in Ref. [44,60]). A similar relation holds for the Z case: Here we did not show explicitly the integration over the invariant mass q 2 of the Z decay particles in the final state, where the Π(q 2 ) dependence enters (see also the Appendix). By taking the ratio of these expressions, the following expression is obtained: In this ratio, any absolute normalization factor appearing in the predicted signal yield, such as the EOT number, cancels out, drastically simplifying the calculation. The same is true for all contributions to the signal efficiency terms η A and η Z that are almost independent from the signal model, as discussed in Sec. IV B. We also observe that this procedure does not depend on the specific value of N up , that can be actually slightly different from the "nominal" setting (2.3 events) due to the inclusion of the non-zero number of expected background events and of the systematic uncertainty factors in the statistical procedure.
To also include the Z resonant production channel in the upper limit calculation, we modified the denominator appearing in the definition of R adding the e + e − → Z → invisible event yield to the total: To solve Eq. 14, we observe that R still contains a residual dependence on g up Z due to the total Z width Γ Z and the additional factor g 2 Z for the e + e − → Z → νν channel (see also the Appendix). To account for this we proceeded by iteration, starting from the ansatz g up Z = (g up Z ) 0 , where (g up Z ) 0 was obtained from the narrow-width approximation ("vanilla" case) or considering the decay to dark sector particles only ("dark" case). At each n−th iteration, we used the value (g up Z ) n−1 to compute Γ Z and obtain (g up Z ) n via Eq. 14. Convergence was observed already after two iterations.

A. Electrons and positrons track-length
We computed the electrons and positrons track-length in the NA64-e ECAL through a Monte Carlo simulation, exploiting the NA64-e GEANT4-based framework [61]. The full NA64-e detector geometry and material composition were implemented in the simulation, including the magnetic field bending the impinging 100 GeV electron beam. Primary electrons were generated just before the upstream tracking stations, with a beam spot size of 1.5 cm and an angular divergence of 0.1 mrad. For all electrons and positrons propagating in the ECAL volume, we sampled the particle energy at each discrete step the trajectory is divided into by GEANT4. We then constructed the corresponding e − and e + energy distributions in the lead and in the plastic scintillator, by assigning to each sampled value a weight given by the step length. The electrons (positrons) track length T − (E) (T + (E)) was finally obtained by normalizing by the total number of simulated events. The obtained result is reported in Fig. 5, displaying the electrons and positrons track length in the lead and in the plastic scintillator. We observe that, due to the ECAL segmentation into equally-sized layers of these materials, the track lengths are almost identical. The factor 2 difference for the high-energy part of the e − distribution is due to the fact that in each layer, including the first, the lead is located in front of the scintillator. Therefore, the energy of the electrons propagating into the first scintillator tile is systematically smaller than that in the first lead layer.
Since, by default, GEANT4 forces a new trajectory step every time a particle crosses the boundary between two regions, we exploited the intrinsic 1.5-mm longitudinal segmentation of the NA64-e ECAL cells to ensure a proper track-length evaluation, without imposing any further subdivision of the particles trajectory. The consistency of the result regarding this choice was checked by repeating the computation of T (E) enforcing a maximum step length of 0.50-mm in the simulation. We observed no significant variations for T + (E). For the T − (E), instead, the two distributions are almost equivalent up to E 80 GeV, while a difference of up to 20% is observed for 80 < E < 99.5 GeV with larger values predicted by the 0.5-mm maximum step-length simulation. For E > 99.5 GeV the difference is even higher, reaching a factor up to 5. The overall normalization of the two distributions is equivalent. We explained this as being related to the aforementioned ECAL geometry. A single 1.5-mm primary electron step in the first lead layer would contribute to the T − (E) sampling with an intermediate E value, whereas if the same step was divided into three 0.5-mm segments, the first would increment the high-energy portion of T − (E). The same effect also applies, in general, for the first ECAL thicknesses, with a reduced intensity. In conclusion, considering that the track-length normalization is not affected by the choice of the stepping size, and that the radiative Z emission has a smooth dependence on the beam energy, in the fol-lowing we will use the T ± (E) result obtained from the nominal GEANT4 simulation.

B. Signal acceptance and detection efficiency
The NA64-e signal efficiency for the A -strahlung channel η A , the Z -strahlung channel η Z ,Rad and the Z resonant production η Z ,Res can be factorized into two different terms. The first, η up , is associated to the response of all detector components installed upstream the NA64e ECAL. This term thus include the tracking efficiency, the efficiency of the SRD cut to reject beam contaminants, as well as the efficiency of the SC counters included in the trigger condition. Global effects such as the overall DAQ efficiency can also be included in this term. The second term η down , instead, is associated to the ECAL, the VETO, and the HCAL detector responses -the main contribution being the 50 GeV ECAL missing energy threshold. By definition, η up is the same for all reaction channels and thus cancels out in the definition of R. Since these channels are characterized by a different signal kinematics, instead, η down has to be computed specifically for each of them.
We evaluated η down for the different reaction channels and parameter models through a GEANT4 simulation of the full NA64-e detector setup using the DMG4 package for events generation [62] and adopting an ad-hoc crosssection biasing mechanism to enhance signal production without distorting the corresponding kinematics, similarly to what was performed in Ref. [60]. The DMG4 package does not offer the possibility to consider off-shell Z production, therefore we used the on-shell approximation Π(q 2 ) → Π(m 2 Z ) in the calculation of η down ; this is justified by the fact that the NA64 detector is sensitive to the missing energy and not to the kinematics of the invisible decay particles. The Monte Carlo event samples were processed through the same NA64-e reconstruction code used for the data analysis, and η down was determined from the fraction of these satisfying all the selection cuts associated to the ECAL, the VETO, and the HCAL, possibly as a function of one or more kinematic observables.
The signal efficiency η down as a function of the Z energy E F is shown in Fig. 6, for the resonant process at m Z = 200 MeV, 250 MeV, and 300 MeV, and for the radiative process at m Z = 3 MeV, 30 MeV, and 300 MeV -these values are representative of the Z mass range explored in this work. We observe that, at fixed E F , all results are compatible with each other within the errors, thus suggesting that kinematic dependence of η down can be effectively taken into account considering its shape as a function of E F . We also checked that the η down A η down Z ,Rad . Therefore, in the following we will use a common expression of η down (E F ) for all reaction channels. The smooth transition observed around E F = 50 GeV is due to the convolution between the 50 GeV threshold on the energy deposited in the ECAL and its finite resolution. Therefore, by including the energy dependence of η down in R, we effectively take into account the modification to the Z line shape due to the detector effects, particularly important in case of resonant production with resonant energy close to the threshold value, i.e. for m Z 2m e E M iss cut 225 MeV.

C. Z events yield
We used the MADDUMP event generator to simulate the Z production in the NA64-e ECAL [63]. MADDUMP is a plugin for the MadGraph5 aMC@NLO program [64,65] developed for fixed thick-target setups that allows to compute the differential yield of Z particles in the lead material of the NA64-e ECAL from the knowledge of electrons and positrons differential track length. For the radiative emission process, we adopted the nuclear form-factor parameterization reported in Ref. [50]. We also explicitly included the factor Π(q 2 ) in the e − e − Z vertex, setting g Z = 1 as justified before. For simplicity, we used an effective polynomial interpolation of the full calculation result presented in Sec. II -to account for the cusp at q 2 = 4m 2 µ , this was implemented separately for the low and high momentum region. Further details are provided in the Appendix.
For a given reaction channel, MADDUMP provides both an unweighted set of N M C Monte Carlo events and the value of the energy-dependent total cross section integrated over the track-length distribution. To include the downstream signal acceptance and detection efficiency, and to account for the ECAL resolution, for each event we computed ε down (E Z ), summed all these values, and normalized the sum to N M C , finally multiplying the integrated cross section by the result. We repeated the calculation independently for the Z -strahlung on the lead nucleus target and for the Z resonant annihilation on atomic electrons. By fixing Π(q 2 ) = 1, we also simulated the radiative dark photon emission on the lead material, necessary to compute the R numerator in Eq. 14.

V. RESULTS
The 90% CL exclusion limits in the m Z vs g Z parameter space obtained from the NA64-e experiment are shown in Fig. 7, for the "vanilla" model (left panel) and for the "dark" one (right panel). In the latter case, to check the effect of changing the Z width, we considered the dark coupling values α D = 0.1 (g D = 1.1) and α D = 0.02 (g D = 0.5), with the fixed mass ratio m Z /m χ = 3. Since, for these values of α D , the missing energy resolution of the ECAL is larger than the Z width, no significant differences are observed between the two cases. Due to tension with perturbative unitarity bound, larger α D values were not considered [25]. The shape of the upper limit curve is associated to the diagram that mostly contributes to the signal yield for a given m Z value. In the region where it is kinematically allowed, roughly defined by E M iss cut < m 2 Z /(2m e ) < E 0 , the resonant emission dominates, resulting in the peculiar "cusp" visible at m Z 250 MeV, while for other m Z values the signal yield is entirely due to the radiative process. The differences between the limit curves in the "vanilla" one and "dark" scenarios are concentrated in the region were the resonant production dominates (∼ 230 MeV −330 MeV). In this mass range the convolution of the Π(q 2 ) function over the Z line shape results in a slight widening of the resonant "cusp" in the "dark" scenario with respect to the "vanilla" one.
We underline that our procedure guarantees that the obtained limit accounts for the effect of all systematic uncertainties that were included in the NA64-e A search analysis (see Ref. [44]). We evaluated the effect of using a global function for η down (E F ) by repeating the calculation of R using separately the expressions for η down Z ,Res , η down Z ,Rad , and η down A ,Rad for the mass value m Z = 250 MeV, where the resonant contribution is dominant. The obtained result is compatible within 1% with the previous one.
We report in the same figure the constraints from other accelerator-based experiments, namely the BaBar search through the visible decay Z → µ + µ − [28] and the Belle-II invisible search result [30]. For the "dark" scenario the BaBar limit, obtained from a search exploiting the Z → µ + µ − decay, does not apply, since for g D g Z the Z decays mostly to the invisible χχ channel. In the same plots, together with the preferred "band" from the muon g−2 anomaly, we also show results obtained by dif-ferent authors through a re-analysis of data reported by neutrino experiments, namely the CCFR result for the trident νN → νN µ + µ − production and the Borexino measurement of solar 7 Be ν e scattering on atomic electrons -however, we point out that these results should be somehow considered cum grano salis, since in both cases not all the experimental details of the original measurement where taken into account in the re-analysis, for example the detector energy resolution; also, the theoretical assumptions for the Borexino limit were questioned in Ref. [36], and a 30% discrepancy was found. The bound obtained from the re-analysis of the COHERENT data, not reported in the figure, is less stringent than those from CCFR and Borexino [37]. We included the sensitivity projection for NA64-µ, a parallel effort of the NA64 collaboration. NA64-µ is a missing-momentum experiment at the CERN M2 beamline, employing the 160 GeV muon beam from SPS to search for the Z via the reaction µN → µN Z and the subsequent invisible Z decay. As discussed in Ref. [42], the sensitivity curve has been computed for an accumulated statistics of 10 11 muonson-target (MOT), for which zero background events are expected. We also report the sensitivity projection for NA64-e for a future high statistics run of 10 13 EOT, assuming the same run conditions of the current e − -beam dataset and considering zero background events. The NA64 collaboration is also investigating the possibility to perform a missing energy experiment with a positron beam, to maximize the signal yield induced by the e + e − channel: we thus show the sensitivity projection for a 10 13 positrons-on-target experiment, again considering zero background events. This result has been obtained following the same procedure used for the electron beam analysis. The track-length and the efficiency were evaluated via GEANT4 and the cross section was numerically integrated with MADDUMP.
The continuous and dashed black curves represent the "thermal target", i.e. the preferred combination of the parameters to explain the observed dark matter relic density, calculated through Eq. 5 by re-scaling the results from Ref. [25].

VI. CONCLUSIONS
The L µ − L τ gauge symmetry extension of the Standard Model provides an elegant explanation to observed "anomalies" between data and SM predictions, such as the muon magnetic moment puzzle. In the "dark" flavour, the corresponding Z gauge boson acts as a portal between SM and a "dark sector", possibly connected with the DM phenomenology. In this work, we presented the exclusion limits for the Z parameters space obtained from the analysis of the existing NA64-e experiment dataset, based on a rigorous treatment of the momentum-dependent coupling between the Z and the first-generation leptons induced by a loop mechanism. These are the first limits set by a direct experimental FIG. 7. The NA64-e exclusion limit for the Lµ − Lτ model, for the "vanilla" (left) and "dark" (right) flavour (red curve). The red (orange) dashed curves represent the sensitivity projections for a future high-statistics NA64-e run with an electron (positron) beam, for a total accumulated charge of 10 13 EOT, while the green dashed curve is the sensitivity projection of NA64-µ [42]. The gray areas are the regions excluded by phenomenological re-analysis of neutrino experiments [34,35], while the blue region is the area excluded by BaBar [28] for the "vanilla" case. Finally, the black curves represent the so-called "thermal target" for the two values of αD = 0.1 and αD = 0.02, i.e. the preferred combination of the parameters to explain the observed dark matter relic density. These have been calculated through Eq. 5 by re-scaling the results from Ref. [25].
search for Z that exclude the region up to the muon g −2 preferred band for m Z 1 MeV, confirming the results already reported by the re-interpretation of neutrino experiments data. Our work demonstrates the potential of the NA64−e experiment, also regarding the complementarity to the future searches with NA64−µ. Future highstatistics NA64-e runs will explore even larger regions in the parameters space, with the positron-beam measurement playing a significant role due to the electronpositron annihilation production mechanism. ACKNOWLEDGMENTS AC and LM warmly thanks Luca Buonocore for very helpful support regarding MADDUMP. AC also thanks Patrick Foldenauer for the useful discussion concerning the connection between the L µ − L τ model and the dark matter phenomenology and Bertrand Echenard for providing the numerical values for the existing BaBar exclusion limits. We warmly thank the reviewer of the manuscript for providing us the full analytical expression for the Π(q 2 ) function. We gratefully acknowledge the support of the CERN management and staff and the technical staffs of the participating institutions for their vital contributions. This result is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, Grant agreement No. 947715 (POKER). This with r µ = m 2 µ /q 2 and r τ = m 2 τ /q 2 .
For simplicity, for the real part of Π(q 2 ) we used an effective parameterization, here denoted by F (q 2 ). Fig. 8 shows a comparison between F (q 2 ) and (Π(q 2 )), nu-merically evaluated: the relative error remains below 5% and (Π(q 2 )) > F (q 2 ) over the considered q 2 range, resulting in conservative limits on g up Z . The imaginary part of Π(q 2 ) can be evaluated analytically, and reads: (Π(q 2 )) =  8. The real part of the Π(q 2 ) function, comparing the numerical result obtained from Eq. 6 with the numerical parameterization FR adopted in this work. For illustration purposes, the arbitrary coupling choice g Z = 2π 2 /e was made.