Confirming $U(1)_{L_\mu-L_{\tau}}$ as a solution for $(g-2)_\mu$ with neutrinos

The recent measurement of the muon anomalous magnetic moment by the Fermilab E989 experiment, when combined with the previous result at BNL, has confirmed the tension with the SM prediction at $4.2\,\sigma$ CL, strengthening the motivation for new physics in the leptonic sector. Among the different particle physics models that could account for such an excess, a gauged $U(1)_{L_\mu-L_{\tau}}$ stands out for its simplicity. In this article, we explore how the combination of data from different future probes can help identify the nature of the new physics behind the muon anomalous magnetic moment. In particular, we contrast $U(1)_{L_\mu-L_{\tau}}$ with an effective $U(1)_{L_\mu}$-type model. We first show that muon fixed target experiments (such as NA64$\mu$) will be able to measure the coupling of the hidden photon to the muon sector in the region compatible with $(g-2)_\mu$, and will have some sensitivity to the hidden photon's mass. We then study how experiments looking for coherent elastic neutrino-nucleus scattering (CE$\nu$NS) at spallation sources will provide crucial additional information on the kinetic mixing of the hidden photon. When combined with NA64$\mu$ results, the exclusion limits (or reconstructed regions) of future CE$\nu$NS detectors will also allow for a better measurement of the mediator mass. Finally, the observation of nuclear recoils from solar neutrinos in dark matter direct detection experiments will provide unique information about the coupling of the hidden photon to the tau sector. The signal expected for $U(1)_{L_\mu-L_{\tau}}$ is larger than for $U(1)_{L_\mu}$ with the same kinetic mixing, and future multi-ton liquid xenon proposals (such as DARWIN) have the potential to confirm the former over the latter.

When combined with the previous Brookhaven determination [26][27][28] this leads to a 4.2 σ observed excess of This exciting result represents a substantial improvement with respect to the previous measurement at BNL [26][27][28], as the experimental uncertainty has been reduced. Since the central value has shifted towards the SM value, the resulting discrepancy falls short of confirming this as evidence for new (BSM) physics, although it certainly leaves the door open for this exciting possibility. While more data from E989 is needed to ultimately judge whether this excess is due to a statistical fluctuation or new physics, complementary search strategies for the possible origin of this excess should be pursued.
During the past two decades, different realisations of particle physics models have been proposed to address this tension, which in general consider extensions in the leptonic sector. Among these, the gauged U (1) Lµ−Lτ stands out for its simplicity, as it is anomaly free without the addition of any extra new fermion fields [29][30][31]. Various realisations of this model have been studied, including also additional neutrino [31][32][33] and dark matter (DM) fields [34][35][36][37]. Recently, it has also been noted that a simultaneous explanation of the Hubble tension can be achieved in these models [38,39]. The purpose of this work is to identify a set of complementary experiments suited for an independent confirmation of U (1) Lµ−Lτ as a solution to the (g − 2) µ anomaly.
The fundamental difficulty in inferring the underlying model of the (g−2) µ excess is that (g−2) µ is only sensitive to the muon coupling of any new hypothetical mediator. Hence, in order to explain the (g − 2) µ excess alone, a simplified mediator model suffices in which the new scalar or vector mediator couples to muons only. Such a new muonphilic singlet mediator will be directly accessible at muon beam experiments. If indeed the only coupling of the new mediator to the SM is to muons alone, it can for example be tested in high-energy muonic Bhabha scattering or µµ → hγ at muon colliders [40][41][42]. However, in order to explain the positive shift in (g − 2) µ with respect to the SM, the new interaction must be (mostly) of scalar-or vector-type rather than pseudo-scalar or axial-vector [43,44]. This necessarily implies that the new mediator has to couple to both left-and right-handed components of the muon with approximately equal strength. Since the fundamental building blocks of the SM in terms of matter fields are the irreducible representations of SU (3) C × SU (2) L × U (1) Y , gauge invariance requires the new mediator at the fundamental level to couple to the left-handed SU (2) L lepton doublet L rather than just the left-handed component of the muon µ L .
In the case of a new vector mediator, the only gaugeinvariant, renormalisable interaction terms required to solve the (g − 2) µ anomaly are given by where L T 2 = (ν µL , µ L ) and g x denotes the coupling of the new vector boson X α . The Lagrangian in Eq. (4) necessarily implies an equal-strength interaction of the new vector boson with the muon-neutrino as with the muon. In this paper, we study the potential of such neutrino interactions for unveiling the underlying physics behind the (g − 2) µ excess, if indeed it is due to a new vector mediator.
The Lagrangian in Eq. (4) could correspond to the interactions of a U (1) Lµ gauge boson with the SM. However, such a U (1) Lµ is anomalous with only the SM field content and new fields need to be added to cancel the anomalies. 2 Sensible UV completions of the U (1) Lµ model quickly run into the highly restrictive flavour-changing constraints [46]. 3 On the other hand, the U (1) Lµ−Lτ model is anomaly-2 Such anomalous leptophilic vector models of gauged lepton number can play an important role in baryogenesis [45,46]. 3 We want to thank Yue Zhang for pointing this out to us. free with only the SM field content and is therefore theoretically very appealing. In this work, we are interested in experimentally confirming U (1) Lµ−Lτ by measuring its three characteristic properties.

P1.
A vector-like coupling to the second generation of leptons. P2. A specific value for the kinetic mixing with the SM photon, namely ∼ g x /70. P3. A vector-like coupling to the third generation of leptons, which is equal and opposite to that of the second.
As explained above, a vector solution to the (g − 2) µ excess only needs to satisfy property P1, so in order to confirm that U (1) Lµ−Lτ is indeed responsible for the (g − 2) µ anomaly, we must verify properties P2 and P3. Since a generic U (1) Lµ -type mediator satisfies property P1, has the freedom to satisfy property P2, but does not satisfy property P3, it is a good model to contrast with U (1) Lµ−Lτ due to its similar, but not identical, physical predictions.
If the new vector boson exists and is light, in particular lighter than the dimuon threshold M X < 2m µ , the only way to directly search for such a boson is via its invisible decay into neutrinos. This can for example be done by a missing momentum search in muon fixed target experiments like M 3 [47] or NA64µ [48,49], where the muonphilic vector is produced via Bremsstrahlung in muonnucleus collisions. Additionally, it can also be looked for in the low-energy e + e − collider experiment BESIII [50]. If on the contrary the new boson is heavy, M X > 2m µ one could directly search for it in the visible decay into a pair of muons at NA64µ [49] as well as at M 3 [47]. Alternatively, the new boson can be produced in kaon decays K → µνX and searched for in a missing energy signature at NA62 [51]. Hence, these searches would be prime candidates to confirm that the (g − 2) µ excess is indeed due to a new muon-philic vector mediator, thus testing P1. However, these missing energy searches are only sensitive to the production cross section of the muon-philic vector boson σ X and its total invisible branching ratio BR X→inv . Therefore, these searches cannot distinguish a muon-coupled boson with decays into light dark sector states from e.g. a U (1) Lµ or U (1) Lµ−Lτ boson with decays into neutrinos.
In order to test P2, we need to get a handle on the kinetic mixing parameter, . To this end, we consider experiments measuring coherent elastic neutrino-nucleus scattering (CEνNS) at spallation sources. This also allows one to distinguish between neutrino-coupled and DM-coupled mediators because the former would induce corrections to the CEνNS rate, while the latter would not. CEνNS was first detected at the COHERENT experiment, which employed a CsI[Na] target [52,53]. Recently, this observation has been replicated using argon nuclei [54]. Since the results are compatible with the SM theoretical prediction, they can be used to constrain new physics in the neutrino sector. Data from the CsI run was used to derive limits for the U (1) Lµ−Lτ model [55] and in Ref. [56] we computed the limits from CENNS-10 LAr results. The combined bounds from CsI and LAr have been presented in Ref. [57]. The current exclusion line does not probe the (g − 2) µ favoured region for U (1) Lµ−Lτ .
In the future, there are plans to build more sensitive multi-ton detectors to further test CEνNS and probe neutrino non-standard interactions [57][58][59]. These include the European Spallation Source (ESS) [60], the Coherent CAPTAIN-Mills (CCM) at Los Alamos National Laboratory [61], and Oak Ridge National Laboratory [62]. In this work, we investigate the sensitivity of all these experiments to U (1) Lµ−Lτ and find that their predicted reach will allow them to test some of the areas of its parameter space which are consistent with the recent (g − 2) µ measurement. In the case of detection, these experiments combined with muon beam experiments provide a powerful strategy for constraining the value of .
It should be noted that all above detectors are only able to probe the second generation lepton couplings, thereby not directly testing P3. In order to confirm a U (1) Lµ−Lτ boson as a solution for (g − 2) µ , one therefore needs sensitivity to the third generation lepton couplings. This can be achieved by considering solar neutrino scattering, as the solar neutrino flux has a ν τ component due to neutrino oscillations. In this context, direct detection experiments are becoming increasingly sensitive to rare events, and they will soon have the potential to probe new physics in solar neutrino scattering. Although their main purpose is to observe dark matter particles, the increased payload, together with the improved cleanliness and analysis techniques, will turn the new generation detectors into multi-purpose instruments. In particular, it is expected that they will be able to observe neutrinos through either their scattering off electrons or through CEνNS. This is, in fact, often regarded as a limitation for DM searches, since CEνNS events will constitute a new background which is difficult to disentangle from potential DM events.
In Ref. [56], we pointed out that direct detection, and in particular experiments based on liquid xenon, could be crucial to probe the U (1) Lµ−Lτ solution to the muon anomalous magnetic moment. The new results of the E989 experiment confirm this, since the 2 σ region in the (g µτ , M A ) parameter space lies partly within the expected sensitivity of future detectors such as DAR-WIN. In this work, we study how well a potential signal from U (1) Lµ−Lτ in agreement with the latest (g − 2) µ measurement could be reconstructed in multi-ton liquid xenon experiments. We find that, since the contribution to nuclear recoils from U (1) Lµ−Lτ and U (1) Lµ differs, direct detection experiments will be crucial to distinguish between these models when combined with results from muon beam and CEνNS experiments.
This article is organised as follows. In Section 2, we introduce two extensions of the SM with a muon-philic vector mediator, namely a minimal gauged U (1) Lµ−Lτ and a simplified U (1) Lµ for comparison. We compute the region of the parameter space compatible with the new (g − 2) µ measurement and define a set of benchmark points for analysis. In Section 3, we study how this region can be probed at muon beam experiments, specifically at NA64µ, arguing that these experiments will be able to constrain the mass of the mediator and the couplings to the muon sector. Next, in Section 4, we turn our attention to experiments looking for CEνNS in spallation sources, reviewing the current experimental situation and the forecast for future detectors, finding that these experiments can constrain the kinetic mixing when combined with muon beam experiments. Finally, in Section 5, we include direct detection probes and study the performance of upcoming and future xenon-based detectors, focusing on the observation of nuclear recoils induced by the CEνNS of solar neutrinos. We show that, when combined with the previous experimental analysis, direct detection experiments could confirm a signal from the U (1) Lµ−Lτ model for mediator masses below ∼ 50 MeV due to their sensitivity to the tau neutrino coupling. Finally, we present our conclusions in Section 6.
2 Muon-philic vector mediators and (g − 2) µ The focus of this paper is to study the physical implications of an explanation of the (g − 2) µ excess in terms of a U (1) Lµ−Lτ gauge boson. For this purpose we consider an extension of the SM by a minimal U (1) Lµ−Lτ gauge model and an effective U (1) Lµ -type model. Throughout this work, we not only test an experiment's ability to reconstruct a signal from a U (1) Lµ−Lτ gauge model, but also test its power to differentiate between the predictions of a U (1) Lµ−Lτ and those of a U (1) Lµ .
Quite generically, the Lagrangian of a U (1) X extensions of the SM can be written in the form where B αβ and X αβ denote the field strength tensors of the SM hypercharge boson B α and new U (1) X gauge boson X α , Y is the kinetic mixing parameter, M X is the mass of the new boson, and g x and j X α are the U (1) X coupling constant and gauge current, respectively. For the two types of models we are considering, the currents are given by where L T i = (ν i , e i ) denotes the lepton doublet of the ith generation and ψ are new heavy fields needed to cancel the anomalies in the case of U (1) Lµ .
It is worth pointing out that effective U (1) Lµ models as an explanation of (g−2) µ as considered here can be subject to very strong constraints from flavour-changing processes like K → πX and B → KX [46]. These constraints are due to enhanced production of the longitudinal mode of the associated gauge boson coupled to an anomalous current [68]. In the low-energy EFT, the heavy fields ψ  [38], Borexino [56], white dwarf cooling [63], COHERENT CsI [55] and LAr [56], neutrino tridents (Charm-II) [64,65] and BaBar [66] and CMS [67] 4µ searches in grey. We show the neutrino trident constraint from CCFR as a grey dashed line, since some backgrounds have not been properly taken into account [51]. The corresponding plot for a simplified U (1)L µ with the same kinetic mixing looks exactly the same with the exception of the BaBar and CMS 4µ limits, which are slightly shifted towards smaller couplings by a factor BR µ−τ A →µµ /BR µ A →µµ ≈ 0.87 (0.71) below (above) the ditau threshold due to the increased branching ratio of the Lµ boson to muons. The green band shows the region favoured at 2σ by the recent confirmation of the (g − 2)µ excess by the E989 experiment [1]. The blue band shows the region favoured by H0 [38]. For illustrative purposes we also show the four benchmark points BP1 -BP4 used for analysis in this work.
(which carry electroweak charge) can be integrated out to produce Wess-Zumino (WZ) terms coupling the new X-boson to the weak bosons, leading to flavour-changing penguin diagrams. The coefficients for the WZ terms are dictated by the underlying UV completion of U (1) Lµ and are, in general, non-vanishing.
Since we are using the U (1) Lµ model as a mere foil for the U (1) Lµ−Lτ model, we are not interested in its possible UV completions, and we rather view it as an effective model for a muon-coupled vector mediator. If indeed experimental results begin to favour U (1) Lµ , this could indicate some finely tuned cancellation of flavour-changing currents. Alternatively, there could be a scenario where the strong flavour constraints are absent. This could be realised with UV-completing heavy fields ψ that are SMchiral with masses that break the electroweak symmetry [68]. The new fields ψ would enter the spectrum at a new physics scale M NP somewhere between the electroweak scale v = 246 GeV and a possible GUT scale f GUT ≈ 10 16 GeV. 4 An extension of the SM by a gauge U (1) Lµ−Lτ is already anomaly-free with only the field content of the SM [29][30][31], but it can even be extended by three right-handed neutrinos to a minimal physically viable model for neutrino masses with correct prediction of the CKM and PMNS matrices [32,[71][72][73][74].
In order to make contact with experiments, we have to canonically normalise the fields in Eq. (5) and change to the mass-basis of the fields. This can be done by diagonalising the kinetic terms of the gauge bosons via a field redefinition and two consecutive rotations outlined in Ref. [63], leading to the interaction terms in the form where e denotes the electromagnetic coupling constant and g Z = e/(sin θ W cos θ W ) with the weak mixing angle θ W . Here j em α and j Z α are the electromagnetic and weak neutral current, and A α , Z α and A α are the mass eigenstates of the SM photon, weak neutral boson and the new gauge boson, which we will refer to as hidden photon. The coupling matrix K is approximately given by [63] with the physical kinetic mixing parameter x = Y cos θ W .
In this article, we will assume that there is no kinetic mixing at tree level, either because a mixing term is forbidden by the underlying UV symmetries or because the tree-level mixing is much smaller than the loop-induced one. At low energies q v EW the mixing of the new Xboson in Eq. (5) is to good approximation only with the SM photon A. In this limit of a low-energy effective theory, the kinetic mixing parameter induced at one-loop can be expressed as where f are the fermions running in loop, Q f , Q x f and m f denote their electromagnetic, U (1) X charges and mass, respectively, and µ is a renormalisation scale. For our models, the kinetic mixing parameters can then be approximated at energies below the muon mass, q 2 m 2 µ , by where the numerical prefactor Log ≈ 10 −2 − 10 −1 for a new physics scale for the fields ψ of M NP = 10 2 −10 16 GeV. So, even if we are agnostic about the true mass scale of the new fields ψ, we can get a rough estimate of the loopinduced kinetic mixing in the case of a UV completion of U (1) Lµ .

Muon anomalous magnetic moment
The muon-philic gauge bosons of U (1) Lµ and U (1) Lµ−Lτ lead to a shift in the anomalous magnetic of the muon 5 , 5 Note that a similar contribution exists for the anomalous magnetic moment of the electron ae, which is, however, doubly where α x = g 2 x /4π, x µ = m µ /M A and Q x µ denotes the charge of the muon under the new gauge group.
Since in both U (1) Lµ and U (1) Lµ−Lτ the charge of the muon is fixed at Q x µ = 1, the observed excess in Eq. (3) translates directly to a band of preferred coupling g x and mass M A values in the 2D parameter as shown in Fig. 1. This defines a clear target parameter space for where to look for this kind of new physics. In this paper, we want to explore the next steps necessary to test the hypothesis that the observed (g − 2) µ excess is indeed due to a U (1) Lµ−Lτ hidden photon by exploring this target region with special emphasis on the additional insights gained by exploiting the neutrino interactions of these bosons.

Benchmark points and analysis strategy
In this work, we want to explore the sensitivities of various different experiments like muon beam, coherent elastic neutrino nucleus scattering (CEνNS) and DM direct detection experiments to a solution of (g − 2) µ in terms of a muon-philic hidden photon A . For this purpose, we define the four benchmark points (BP) of Table 1 in the region of parameter space favoured by (g − 2) µ in a U (1) Lµ−Lτ model. The benchmark points are also illustrated in Fig. 1.
Throughout this paper, we study the sensitivities of the different experiments to a hypothetical signal coming suppressed by kinetic mixing. Since we are only considering loop-induced kinetic mixing, this corresponds effectively to a three-loop process. In this context, it is interesting to note that the latest experimental determination of (g − 2)e shows a mild deficit compared to the SM prediction [75]. Due to the positive sign of the shift from vectorial couplings, the corresponding contribution via a U (1)L µ or U (1)L µ−Lτ boson slightly worsens this tension.
from a U (1) Lµ−Lτ boson with coupling and mass as defined by the benchmark points. For all experiments, we therefore generate mock data sets for the four different U (1) Lµ−Lτ benchmark points and try to reconstruct the mass M A and coupling g x under both the hypothesis that the signal is due to a U (1) Lµ−Lτ or due to a U (1) Lµ hidden photon. In the remainder of this paper, we will refer collectively to the coupling and mixing of U (1) Lµ and U (1) Lµ−Lτ with g x and x , and label them by an index µ or µτ when we need to distinguish between these groups.

Muon beam experiments
With the recent measurement of the E989 experiment confirming the (g − 2) µ excess, independent complementary measurements will be of paramount importance for pinning down its true nature. If the excess is indeed due to a new light muon-philic vector boson, a missing energy search at muon beam dump experiments will be crucial in confirming this. Both the planned CERN NA64µ [78,79] and Fermilab M 3 [47] experiments are prime candidates to conduct such a search. The NA64µ experiment is planned to perform a first pilot run in 2021 [80] and conduct its full Phase-I run in 2021-2023 [79] with a total of 10 11 muons on target (MOT). This will therefore be the first experiment to be sensitive to the full parameter space of a muon-philic boson allowed by (g − 2) µ . We therefore constrain our analysis of missing energy searches at muon beam experiments to the case of NA64µ, but note that our results should be easily generalised to the case of the Fermilab M 3 experiment.
Complementary to these two muon fixed target experiments, is the search for invisible decays of a muonphilic mediator A in rare kaon decay experiments like NA62 [51]. There, the A could be produced from the final state muon in the kaon decay K → µν. The A can then decay invisibly, leading again to a missing energy signature. While in principle this experiment is also capable of probing the (g − 2) µ region, NA62 will need its total planned number of 10 13 collected kaons. More importantly, however, NA62 is not a background free experiment, and in order to be sensitive to U (1) Lµ−Lτ as a solution to (g − 2) µ , the experimental systematics have to be further reduced [51]. We leave a detailed analysis of NA62 for future work and focus on NA64µ in this paper. 3.1 Signals of (g − 2) µ at NA64µ In the NA64µ experiment, a beam of muons with E 0 ∼ 160 GeV is dumped onto a lead target with a thickness of L T = 40 X 0 ∼ 20 cm. The energy and momentum of the scattered particles are measured in the fiducial volume of the detector with L D = 5 m, consisting of the active target ECAL, a magnetic spectrometer, tracker and HCAL [79]. The production of a light muon-philic vector boson in the NA64µ target proceeds via the Bremsstrahlung process µ + Z → µ + Z + A illustrated in Fig. 3. Since the virtuality of the photon exchanged between the muon and the nucleus is quite small [81], the 2 → 3 production cross section is related to real photon scattering via the Weizsäcker-Williams approximation [82,83], where α is the fine structure constant, A /E 2 0 is the boost factor of the hidden photon, x = E A /E 0 is the energy fraction carried away by the hidden photon, and χ denotes the effective photon flux sourced by the nucleus. The photon flux is given in terms of the electric form factor G 2 (t) of the nucleus as [82,84] where t = −q 2 and one can approximate t min = (M 2 A /2E 0 ) 2 and t max = M 2 A . The electric form factor has an elastic and an inelastic contribution, which are given e.g. in Ref. [81]. Taking the expression for real photon scattering from Ref. [48] and integrating over the hidden photon emission angle θ, the 2 → 3 differential cross section in the Weizsäcker-Williams approximation is where and is the virtuality of the intermediate muon.
The total number of expected hidden photon events at NA64µ can then be expressed as [85] where y is the penetration depth of the muon in the target, n atom denotes the number density of nuclei in the target material, and x min is a cut imposed on the minimum energy fraction carried away by the hidden photon in order to suppress the background. P (z) denotes the decay probability density of the hidden photon and can be written as where is the decay length. However, as long as the hidden photon is lighter than the dimuon threshold, M A < 2m µ , it will only decay invisibly (except for a very suppressed decay into e + e − pairs) and BR(A → inv) ≈ 1. In this regime, we can effectively set z min = 0 and z max = ∞ so that the integration over P (z) yields one. Following Ref. [85], we note that, at the relevant energies of a few GeV to ∼ 100 GeV, the muon is a minimum ionizing particle with a rather flat stopping power dE/dy so that we can assume it to be constant [86]. This can be used to rewrite Eq. (21), via in terms of the muon energy E µ as Integrating dE/dy over the length of the target L T ≈ 40 X 0 = 20 cm, one finds that the muon hardly loses any energy, and therefore the expression can be further simplified to [85,87] where we have expressed the number density n atom = ρN A /A in terms of Avogadro's number N A , the density ρ, and the mass number A of the target material.
In the case of M A < 2m µ for both U (1) Lµ and U (1) Lµ−Lτ , the total number of invisible A decays is given by Eq. (26). However, in the regime where M A ≥ 2m µ , we have to be more careful since the muon-philic A can also decay visibly into muons. The invisible branching fraction is now altered, but also dimuon decays that occur outside the fiducial detector volume will be invisible. In this case, the expression for the total number of invisible events is given by where L det = L T + L D is the length of the target plus detector. In the relevant coupling regime of g x ∼ O(10 −4 ), the decay of the A is prompt at NA64µ for M A ≥ 2m µ so that almost all of them already decay within the target. The missing energy signature at NA64µ is defined by E tot = E ECAL + E HCAL 12 GeV and a single muon carrying roughly E µ E 0 /2 of the initial energy E 0 . Detailed studies of the expected backgrounds have shown that imposing a maximum muon energy of E µ 100 GeV makes this search essentially background free [48,79]. In the NA64µ proposal, the momentum resolution is given as σ p ∼ 1.3 GeV and σ p ∼ 2.3 GeV for a measurement of the incoming muon momentum at the beam muon station and magnetic spectrometer, respectively [79]. We assume that the resolution of the measurement of the outgoing muon momentum will be of similar order. Hence, we assume a minimum bin width of 5 GeV for the expected signal counts as a function of the reconstructed muon energy. In order to guarantee that we have roughly O(3) events per bin for our benchmark points, we choose a more  Table 1. The grey shaded areas show the experimental energy cut applied in the measurement in order to remove any backgrounds.
conservative binning of 10 GeV steps for our analysis. Following Ref. [48], we also assume a signal window of 10 GeV to 100 GeV in the muon momentum. In Fig. 4, we show the projected muon spectra produced by a muon-philic vector boson in the Phase-I run at NA64µ for the benchmark points defined in Table 1. These results are obtained with an assumed average signal reconstruction efficiency of = 0.3. Since, to our knowledge, no energy-dependent efficiency curves have been published, this is an approximation of realistic efficiencies, which lie in the range of = 0.1 − 0.5 for masses of 1 MeV -1 GeV [79]. It can be seen that all four benchmark points produce a clearly visible signal with O(1 − 10) events in each bin within the signal region.

Parameter estimation after future detection
Finally, with these projected signal spectra for our benchmark points at hand, we are interested in the question of how well the mass and coupling of such a muon-philic hidden photon can be reconstructed after a potential discovery at NA64µ. In order to answer this question, we perform a maximum likelihood parameter estimation. We build our binned likelihood function based on Poisson distributed events, where N is the total number of bins, n i are the observed number of events in bin i, and µ i = s i + b i is the combined number of expected signal s i and background events b i for the model parameters θ = (g x , M A ). Since the NA64µ signal region is chosen such that the measurement is background free, we can take b i = 0.
Finding the set of parameters θ 0 that maximise the likelihood in Eq. (28) allows us to construct the log likelihood ratio for parameter estimation according to [88], Since under random fluctuations in the data n i the function −2 ln λ(θ) asymptotically follows a χ 2 distribution with the number of degrees of freedom n equal to the number of estimated parameters θ, we can reject parameters θ at the level of 1 − p via [89], where F χ 2 denotes the cumulative χ 2 distribution with n degrees of freedom. Since we are estimating two parameters, the mass M A and the coupling g x of the new muonphilic boson, we can exclude parameters at the 68% CL and 95% CL via −2 ln λ(θ) > 2.30 and −2 ln λ(θ) > 6.18, respectively.
Using this method, we perform a 2D global parameter scan using the expected signal spectra of a U (1) Lµ−Lτ hidden photon at the benchmark points of Table 1 as mock data. The resulting contours in the (g x , M A ) parameter space from both reconstructing these signals within U (1) Lµ−Lτ and U (1) Lµ are shown in Fig. 5.
Let us note that NA64µ has very good sensitivity to mid-range hidden photon masses, M A ≈ 100 MeV, in the (g − 2) µ region. This corresponds to the situation at BP4, where the 68% and 95% CL contours both close around a rather small region enclosing the true parameter point. For the lower mass BP1 -BP3, NA64µ can reconstruct the coupling g x quite well while it can only give an upper limit on the mass. This behaviour can be understood qualitatively from the spectral shapes of the signal. The spectral shape of the scattered muon is entirely fixed by the mass of the new mediator, M A . The coupling g x can thus be viewed as a scaling parameter that determines the overall normalisation of the spectrum. As can be seen in Fig. 4, BP4 leads to a falling spectrum (with increasing momentum) in the signal window. This flattens out slightly between ∼ 60 − 90 GeV before it falls off again more steeply towards higher momenta. This characteristic flattening is a visible feature of the spectrum within the signal window, leading to a good reconstruction of the parameters. However, BP1 -BP3 lead to a rising spectrum with their maxima lying outside the signal window at high momenta. Within the signal window, their spectral shapes are relatively similar with the characteristic peaks lying outside the signal window. This behaviour persists for smaller mass bosons, making their spectra very hard to distinguish from one another.
In general, it can be seen that the signal can be equally well reconstructed within U (1) Lµ−Lτ and U (1) Lµ , and the resulting contours exactly coincide for M A < 2 m µ (i.e. for  all our BPs). For masses M A ≥ 2 m µ , the contours for U (1) Lµ and U (1) Lµ−Lτ would deviate due to their invisible branching fraction of BR inv ≈ 0.33 and BR inv ≈ 0.5 for U (1) Lµ and U (1) Lµ−Lτ above the dimuon threshold, respectively. The reason why the contours exactly coincide for M A < 2 m µ is due to the fact that in this regime the spectra produced by a U (1) Lµ or U (1) Lµ−Lτ A are exactly identical since the production cross section only depends on the muon coupling and both A have an invisible branching fraction of BR inv ≈ 1. Note that, while NA64µ in general reconstructs the coupling g x quite well, it is not at all sensitive to the kinetic mixing parameter x (as long as x g x ). This makes it in general impossible to distinguish a potential U (1) Lµ signal from U (1) Lµ−Lτ by using NA64µ data alone. Hence, complementary probes of such a potential signal are needed to discriminate the two models.

CEνNS at spallation sources
New physics in the neutrino sector induces contributions to the coherent elastic scattering of neutrinos off nuclei. After the first measurement of this rare SM phenomenon by the COHERENT collaboration at a spallation source, the potential of this type of experiment to probe nonstandard neutrino interactions has been thoroughly studied [55,56,58,59,[90][91][92][93][94][95][96][97][98][99][100]. These experiments benefit from the high neutrino flux that can be achieved at spallation sources and the scalability of their detectors. In particular, a number of ton and multi-ton liquid argon detectors have been proposed that will considerably improve the sensitivity to light mediators. At spallation sources, protons collide on a nuclear target to produce pions. These decay at rest to produce a monochromatic 29.8 MeV muon neutrino flux (π + → µ + ν µ ), while the subsequent decay of the anti-muon produces a delayed beam of electron neutrinos and muon antineutrinos (µ + → e + ν eνµ ). The corresponding fluxes (see e.g., Refs. [53, 101]) read The normalisation factor η = rN POT /4πL 2 takes into account the number of neutrinos of each type produced from each proton on target (POT), r, and the baseline, L. The number of expected CEνNS events can then be expressed as where (E R ) is the energy-dependent efficiency, and the minimum and maximum neutrino energies are E min ν ≈ M N E R /2 and E max ν = m µ /2. The number of target nuclei is N targ = M tot /M N , and M N the atomic mass (we will concentrate on LAr detectors, for which we will assume 100% isotopic abundance of 40 Ar). The CEνNS differential cross section reads [56], where the coherence factor for the effective neutrinonucleus coupling via the SM Z-boson is defined as Q νN = N − (1 − 4 sin 2 θ W ) Z, and F (E R ) is the Helm nuclear form factor [102,103].
Since there is no tree-level coupling between the A and first generation leptons in either U (1) Lµ−Lτ or U (1) Lµ models, and because the ν τ flux is negligible (see e.g. Ref. [58] for upper limits from oscillation parameters), the new physics contribution to N CEνNS only probes the coupling between A and the muon sector (α = µ in Eq. (33) and Eq. (32)). Also, as the muon neutrino charge, Q x νµ , is positive and x < 0, the new physics contribution from both U (1) Lµ−Lτ and U (1) Lµ results in a negative interference term. For the relevant part of the parameter space, this leads to a reduction in the number of expected events with respect to the SM prediction 6 , which makes these scenarios more difficult to probe.
In this section, we use data from the benchmark points defined in Table 1 to obtain a reconstruction in the (g x , M A ) plane using future CEνNS detectors. Since the number of events is a function of the product g x x , the predictions from U (1) Lµ−Lτ and U (1) Lµ are, in principle, indistinguishable from each other. The reconstructed value of g x can be interpreted as different values of g µτ and g µ (obtained using the corresponding relation for x from Eq. (11) and Eq. (12), respectively). However, one can expect that the combination with other types of experimental searches would help to discriminate these possibilities, as we will comment below.
We have considered four detector configurations, based on planned CEνNS experiments, with parameters as detailed in Table 2, inspired by the analysis of Ref. [58]. First, we include a 610 kg (fiducial mass) extension of the current CENNS-10 LAr detector [62] at the SNS. In order to simulate the number of expected events, we have considered the same quenching factor to relate the nuclear and electron equivalent energies ( ). The efficiency is expected to be similar to that of its predecessor, which dropped to 50% at an energy of approximately 4 keV ee (equivalent to approximately 20 keV nr ). We have approximated it by = 0.5 (1 + tanh(E R − 4 keV ee )) and obtained good agreement with the spectrum predicted in Ref. [58]. Regarding the systematic uncertainty, we have fixed it at 8.5%, as in CENNS-10. Based on the projected characteristics of CCM [61], we have considered a fiducial mass of 7 tons of LAr. CCM is planned to run for a total of 2.5 years at Lujan in a near position (L = 20 m) and a far position (L = 40 m) configuration. We have only used data from the near position run (as the expected number of CEνNS events is approximately four times larger), assuming a total operation time of 1 yr in this configuration. Regarding the ESS [60], we proceed as in Ref. [58] and consider two different setups: a small (10 kg) but extremely sensitive (E th = 0.1 keV) detector, and a large one (1 ton) with the same threshold energy as CCM and CENNS. For both configurations, the baseline is L = 20 m and we assume 1 yr of operation.
For both CCM and ESS, we consider 100% efficiency and we do not include any quenching factor. Regarding the threshold, we will assume a sharp cut at E th = 20 keV nr (this value is similar to the energy at which the efficiency of CENNS-10 LAr drops to 50%).
An important background for this kind of search is due to beam related neutrons (BRN), which for CENNS-10 represented approximately 10% of the signal events [62]. Since we ignore how CCM and ESS would perform in this respect, we have assumed N bkg = 10%N SM . We  Table 1 using CEνNS data from the experiments of Table 2 and fixing x = −gx/70. The solid (dashed) contours correspond to the 68% CL (95% CL) and the star represents the position of each benchmark point. The dot-dashed lines denote the 90% CL upper limit on the coupling gx in the case no reconstruction is possible.
will consider four energy bins across the energy window of 20 − 100 keV. Having such a modest measurement of the energy spectrum allows for some limited reconstruction of the mass of the mediator in the event of an observation, but most importantly, it helps reducing the effect of the normalisation systematic error.
We have used a χ 2 test to compare the observed number of events in the i-th energy bin, N i obs , with the theoretical prediction for a given point in the parameter space, N i th (g x , M A ), The statistical uncertainty is defined as where N i bkg = N i SM /10 (and is therefore the same for all benchmark points). There is a systematic uncertainty in the total normalisation, represented by a nuisance parameter, a, which we take as 5%. We consider this is a realistic goal, based on the 8.5% systematic error achieved in CENNS-10 LAr.
The energy spectra for the four benchmark points of Table 1 are represented in Fig. 6, where the SM prediction is shown as a black solid line for comparison. As mentioned above, all BPs show a deficit of events with respect to the SM rate. The difference is much more prominent at low energies, which is an excellent motivation to reduce the experimental threshold in this kind of experiment. ference term dominates). This means that only the lowmass window of the area compatible with (g − 2) µ can be probed: of all the benchmark points in Table 1, BP1 and BP2 would be observed in the future ESS. For BP3 and BP4, we obtain upper bounds on the coupling as a function of the mediator mass.
In particular, the upcoming CENSS610 and CCM experiments would be unable to disentangle the SM CEνNS flux from a new physics signal, resulting in upper bounds on the parameter space that reach couplings of the order of g x ∼ 7 × 10 −4 (for CCM, this is in agreement with the results from Ref. [57]). The predictions for ESS are more optimistic. The case of ESS10 perfectly illustrates the great benefit of lowering the experimental threshold: a small but extremely sensitive configuration would be able to observe the nuclear recoil spectrum where the new-physics con-tribution from light mediators is maximal (see the lower left panel in Fig. 6). Not only does this allow the new physics prediction to be disentangled from that of the SM (thereby producing closed contours in g x ), but also some sensitivity to the mediator mass is obtained (and upper bounds are found for M A ). The bump at low mediator masses in the ESS10 lines corresponds to the crossover of the contribution from the interference term and the pure BSM one (only observable at low energies). The 1 ton configuration of ESS would perform slightly better at larger masses, and could confirm observation for BP1 and BP2 (providing closed contours in the coupling g x but a poor reconstruction of the mediator mass).
It should be emphasised that the results of Fig. 7 can be interpreted in both the U (1) Lµ−Lτ or U (1) Lµ model. Experiments looking for CEνNS at spallation sources, on their own, would be unable to measure the kinetic mixing parameter x .

Combination with NA64µ
The combination of results from NA64µ and CEνNS experiments at spallation sources can potentially shed some light on the value of the kinetic mixing parameter. As explained in Section 3, NA64µ is only sensitive to g x , whereas CEνNS probes g x x . Thus, the combination of results from both sources could be used to infer the value of x . In this section we illustrate this complementarity, showing that tension in the reconstructed regions arises if the wrong assumption is made for x .
First, Fig. 8 shows the combined results from future spallation sources and NA64µ, assuming x = −g x /70. This coincides with the prediction in a U (1) Lµ−Lτ model and with the assumption that we made in generating the benchmark points. The evidence strengthens for BP1 and BP2, since these benchmark points can be probed by both types of techniques. Even if spallation source experiments do not observe new physics (as is the case of BP3 and BP4) the resulting bounds can narrow down the region compatible with NA64µ, providing a better measurement of the mediator mass.
In contrast, Fig. 9 assumes x = −g x /10 in the reconstruction, a value that could come from a particular realisation of a U (1) Lµ model. As expected, the reconstructed areas show increasing tension. For example, for BP1, there is a small overlap of the 95% CL contours, which shrinks for BP2. In contrast, for BP3 and BP4 almost the whole area compatible with NA64µ is excluded, indicating a wrong assumption for x . Measuring x does not necessarily imply discriminating U (1) Lµ−Lτ and U (1) Lµ . In particular, if x is shown to be inconsistent with −g x /70, then the minimal U (1) Lµ−Lτ model would be ruled out as an explanation for (g − 2) µ . However, if x is found to be compatible with −g x /70, this could correspond to either U (1) Lµ−Lτ or U (1) Lµ . The only way to tell these constructions apart would be to explore the couplings to the tau sector. As we will argue in the next section, this could be done by observing solar neutrinos in direct detection experiments.

Direct detection experiments
The expected differential rate of solar neutrino-electron or CEνNS scattering events at direct dark matter detection experiments can be written as is the minimum neutrino energy to produce a nuclear recoil or electronic recoil of energy E R , m T is the mass of the target, and n T is the total number of targets (electrons or nuclei) per unit mass. Regarding the solar neutrino fluxes, dφ νe /dE ν , we will consider a high metallicity scenario, which increases the flux of neutrinos in the pp-chain by up to 10%. For the case of 8 B neutrinos, this corresponds to a total flux of 5.46 × 10 −6 cm −2 s −1 [104]. The energy-dependent oscillation probability P (ν e → ν α ) is taken from Fig. 7 of Ref. [56].
The differential cross section of the elastic scattering of neutrinos off nuclei is given in Eq. (33), and the corresponding expression for scattering off electrons for α = µ, τ reads where g e L = sin 2 θ W − 1/2 and g e R = sin 2 θ W are the couplings of the Z-boson to the electron. The SM expression for α = e is obtained by setting x = 0 and replacing g e L → 1 + g e L in the above expression. For more details, see Ref. [56]. The first term in Eq. (36) corresponds to the pure SM, the second to the interference, and the last one to the pure BSM contribution. The sign of the interference term critically depends on the charge Q x να of the scattered neutrino ν α , which is positive for muon and negative for tau neutrinos. For electron recoils, the pure BSM term dominates the scattering cross-section at couplings relevant for explaining (g − 2) µ , and therefore we expect an increase in the number of events for both U (1) Lµ−Lτ and U (1) Lµ .
When calculating the expected number of scattering events in direct detection experiments, we include detector-specific effects, such as detector efficiencies and resolutions. Efficiencies (E R ) are folded into Eq. (35) before the theoretical differential rate is convoluted with a Gaussian resolution function with energy-dependent width, σ(E R ). This results in an expected observed count of where ε is the exposure of the experiment, and E th is its energy threshold. An important difference with respect to the physics at NA64µ and CEνNS experiments at spallation centres is that there is now an incoming flux of tau neutrinos. In fact, given the neutrino oscillation probabilities of Ref. [56], the flux of ν τ is larger than that of ν µ at the energies relevant for the 8 B and hep solar fluxes. This is extremely relevant for discriminating new physics models. In particular, the U (1) Lµ−Lτ contribution for both NR and ER is positive with respect to the SM. In contrast, the U (1) Lµ model with µ = −g µ /70 leads to a negative contribution to NR (from the dominant interference term) and a positive one for ER (from the pure BSM term). Therefore, an observation of NR in direct detection experiments would unequivocally point to either a U (1) Lµ−Lτ or U (1) Lµ model. This is by no means an easy task: although the flux of solar neutrinos is much larger than that of atmospheric ones, their energies only reach approximately 15 MeV (for 8 B). These lead to ∼keV nuclear recoils, right at the threshold of current liquid xenon detectors. On the other hand, although solid-state detectors benefit from a much lower threshold, their planned target size makes them unlikely to observe enough events from solar neutrinos in the near future. Thus, in our analysis, we will concentrate on planned LXe detectors and will base our experimental configurations on the upcoming multi-ton LZ [105] and XENONnT [106] detectors and on the projected DARWIN observatory [107]. The proposed exposures are 15.34, 20 and 200 ton · yr for LZ, XENONnT and DARWIN, respectively.
For LZ, we have taken the efficiency functions given in [108], the resolution fit given by LUX in [109], and ER backgrounds from [110]. For XENONnT, the NR and ER efficiency functions have been taken from [111] and [112], respectively, the resolution function has been read from [112], and the ER backgrounds have been taken from [106]. Finally, for DARWIN, we have assumed the same efficiency and resolution functions as for XENONnT, used the lower ER background predictions given by [113], and assumed a flat background rate for all background components except for the double-β decay of 136 Xe. For all experiments, we have taken the nominal energy thresholds to be at the energies where their (event-specific) efficiencies reach 50%. This corresponds to LZ thresholds of 3.8 keV nr and 1.5 keV ee , and XENONnT and DARWIN thresholds of 5.7 keV nr and 1.5 keV ee .
We have performed three different types of analyses: an NR-only analysis (based on the discrimination of nuclear and electron recoils to reduce background), an ERonly analysis, and a combined NR + ER analysis. For the NR-only analysis, we assume a 50% acceptance cut above the energy threshold after rejecting ER events and treat this as a background-free analysis. For the ER-only analysis, we assume 99.5% ER/NR discrimination and include experimentally-relevant ER backgrounds. For the NR + ER analysis, we combine both types of events, interpreting all events as ER events, with the ER background included but without cuts. In each analysis, the appropriate NR/ER efficiency functions are folded into the differential rate, with the energy resolution taken at the appropriate energy scale.
The left panel of Fig. 10 shows the predicted differential rate for nuclear recoils in a generic xenon detector. The NR signal from solar neutrinos from the 8 B flux falls abruptly at approximately 1 − 2 keV. This makes their observation very challenging, forcing experiments to achieve a very low energy threshold. For the nominal value of the threshold of the planned LZ and XENONnT detectors (of the order of 3 keV), only a handful of these neutrinos can actually be observed (mostly due to the effect of the resolution near the threshold). In order to exploit this signal one must therefore achieve lower experimental thresholds. The figure also illustrates how the U (1) Lµ−Lτ model predicts an excess of events with respect to the SM, whereas U (1) Lµ would lead to a reduction of the rate. Given the importance of nuclear recoils in discriminating between U (1) Lµ and U (1) Lµ−Lτ , we take the liberty to vary the threshold, exploring their effect with regards to discovery potential.
The right panel of Fig. 10 corresponds to the expected electron recoil rate for BP1, compared to the SM prediction. For comparison, we also show the expected background, using the projection for DARWIN as guideline. The electron background is dominated by 136 Xe double-β decay at high energies. For the energy range considered in the analysis, one must take into account the fact that electrons are originally bound to the atom. This results in a series of steps in the spectrum, which in a first approximation, correspond to the different ionisation energies and effectively reduces the number of electrons available at low energies. A more careful treatment in terms of a relativistic random phase approximation (RRPA) leads to a further suppression in the 0.25 − 30 keV ee window according to Fig. 2 of Ref. [114], which we have implemented using an energy-dependent scaling 7 . Note that, below 0.25 keV ee , we have reverted back to the step-function approximation, which acts as an upper bound to the expected rate in the absence of numerical solutions to the RRPA at low energies [114].
We have only shown the predictions for BP1, as that is the most optimistic case. The differential rates for BP2  are a bit smaller and for BP3 and BP4 they lie very close to the SM, being much more difficult to separate. In order to determine the potential of direct detection experiments to disentangle U (1) Lµ−Lτ and U (1) Lµ , we have computed the number of expected events using the experimental configurations of future liquid xenon detectors for the U (1) Lµ−Lτ model in the four benchmark points of Table 1. We find that for the upcoming LZ and XENONnT, the number of counts in both the NR and ER channels will be too small to be able to probe any of the benchmark points at the 5 σ level.
We have then considered the exposure and energy threshold as free parameters and studied the minimal conditions under which a 5σ observation of U (1) Lµ−Lτ could be claimed at each benchmark point. To compute the significance, we construct the (one-tail) p-value, where b is the number of background events, and n obs is the observed number of events (in our case the theoretical expectation for each BP). A one-tail statistic is used as we expect the number of observed events to be only in excess of the background-only result. For a discovery-level significance of 5σ, we require a p-value of 2.87 × 10 −7 .
In our analysis, we solve for the required thresholdexposure pairs needed to produce a discovery-level measurement of U (1) Lµ−Lτ for each of our benchmark points. We take the background-only hypothesis to contain all of the counts from SM CEνNS and experimental background events. We repeat the analysis individually for each type of search: NR, ER and NR+ER. When lowering the experimental threshold, we assume each experiment's respective efficiencies would improve in a way akin to an extension of their original efficiency functions in log-space, such that their new, lower thresholds would be where the now extended efficiency function reaches 50%. In all cases, we take the maximum of the energy window to reside at 30 keV ee (∼ 13 keV nr ), above which the double-β decay of 136 Xe is expected to dominate over the solar neutrino signal [105].
When lowering the thresholds, we have been careful not to allow the experimental resolutions to increase to arbitrarily large values at low energies, making a reconstruction of low-energy thresholds unrealistic. We have therefore capped each resolution function at its value at the original 50% threshold, shifting this capped resolution linearly when lowering the thresholds. To ensure a valid conversion of the resolution functions into NR energy scales, we have taken 0.7 keV nr (∼ 0.1keV ee ) as a minimum threshold, corresponding to the lowest energy at which the Lindhard model has been experimentally verified for LXe detectors [115]. Fig. 11 shows the resulting values of the exposure and threshold energy for which a 5σ discovery of U (1) Lµ−Lτ could be claimed. The area above the dashed blue line would be accessible via NR, the area above the green dashed line corresponds to ER and the region above the red line can be observed via NR+ER. This figure illustrates that the optimal strategy to observe NR is to decrease the threshold, whereas to detect ER one simply needs a bigger target. In the case of BP1 (the only bench-   Table 1 assuming a U (1)L µ−Lτ model. We use direct detection data from a hypothetical low-threshold search (E th = 1 keV) in the future DARWIN detector. For BP1, the black solid (dashed) contours correspond to the 68% CL (95% CL) contours. For the rest of the BPs, the black dot-dashed lines represent the 90% CL upper limit on the coupling gx. The rest of the lines follow the same convention of Fig. 8. The black star marks the benchmark point.
help reconstruct the parameters of the U (1) Lµ−Lτ model. Fig. 12 shows the reconstruction of parameters in the (g x , M A ) plane when direct detection data is added to the results from searches at NA64µ and spallation sources, when these are interpreted in terms of a U (1) Lµ−Lτ model. In line with our analysis of Fig. 11, the 5 σ discovery of BP1 or BP2, when combined with data from NA64µ and spallation sources, would considerably narrow down the area compatible with (g − 2) µ , reducing the uncertainty on both the coupling and the dark photon mass. Notice that, if we were to interpret the direct detection results in terms of a U (1) Lµ scenario, the whole parameter space for BP1 and BP2 would be ruled out 8 . For BP3 and BP4, direct detection will not have enough sensitivity to claim detection. Nevertheless, the bounds derived from direct searches, would further constrain the parameter space. Similar to CEνNS bounds, direct detection will be more effective in constraining mediator masses below approximately 50 MeV in the (g −2) µ region. These bounds will be similar for both U (1) Lµ−Lτ and U (1) Lµ scenarios and no extra information would be gained in that respect.

Conclusions
In this article, we have explored the complementarity of different experimental probes of muon-philic solutions to the 4.2 σ deviation between the observed value of the muon anomalous magnetic moment at the E989 experiment and the SM theoretical prediction. In particular, we have focused on the light vector mediator that arises from the anomaly-free U (1) Lµ−Lτ model. We have laid out a strategy of how to combine muon beam experiments, CEνNS experiments at spallation sources and DM direct detection experiments to confirm whether the observed (g − 2) µ excess is indeed due to a U (1) Lµ−Lτ boson. We have done so by contrasting our findings for U (1) Lµ−Lτ to those of a phenomenological U (1) Lµ model, which can experimentally mimic many but not all of the properties of an U (1) Lµ−Lτ .
We can summarise our findings as follows: • NA64µ will be the first experiment to be sensitive to the (g − 2) µ region in U (1) Lµ and U (1) Lµ−Lτ . For mediator masses of the order of M A ∼ 100 MeV, it would allow for an excellent reconstruction of both the coupling g x and the mass M A . For lighter hidden photons (M A 50 MeV) it can still give a good reconstruction of the coupling while only providing an upper bound on the mediator mass. In general, it is, however, insensitive to the kinetic mixing x . Furthermore, with data from NA64µ alone the models U (1) Lµ−Lτ and U (1) Lµ cannot be discriminated (nor can they be distinguished from a muon-coupled mediator that can decay into light DM states).
• Future experiments looking for CEνNS at spallation sources, such as the planned CENNS610 and CCM, will set constraints on the low-mediator mass (M A 50 MeV) region of the (g − 2) µ solution. The predicted CEνNS rate for the U (1) Lµ−Lτ and U (1) Lµ models is smaller than in the SM, making this observation more challenging. Large detectors near powerful sources such as the projected experiment at the ESS or smaller devices with an extremely low-threshold (we consider a E th = 0.1 keV version of ESS with just 10 kg) could be able to reconstruct the coupling and set upper bounds on the mediator mass. However, since these experiments only test the coupling of the hidden photon with the muon sector through the combination g x x , they will be unable to discriminate the U (1) Lµ−Lτ and U (1) Lµ models on their own.
• From the combination of data from NA64µ (which measures g x ) and CEνNS experiments (which measure g x x ), one can infer the value of the kinetic mixing and also improve the reconstruction of the mediator mass. However, this would be insufficient to distinguish the U (1) Lµ−Lτ and U (1) Lµ model as, in principle, the kinetic mixing could be the same ( x = −g x /70).
• Direct detection experiments will provide a unique opportunity to test the couplings to the tau sector via scattering off solar ν τ 's. In particular, because of the contribution of tau neutrinos, the expected nuclear recoil rate from CEνNS in U (1) Lµ−Lτ is larger than the SM prediction, whereas for the U (1) Lµ with µ = −g µ /70 a decrease is expected. This allows for discrimination of both solutions to (g − 2) µ in the values of the kinetic mixing to which NA64µ and spallation experiments are most insensitive. We have estimated the experimental requirements to probe the low-mass window of the (g − 2) µ region. We conclude that a liquid xenon experiment with a total exposure of 100 ton yr and a threshold of E th = 1 keV would suffice to probe mediator masses below 25 MeV. These are similar characteristics to those of the planned DAR-WIN detector. The mass and coupling reconstruction in that area would be complementary to that of NA64µ and CEνNS experiments.
In the event of no observation, DARWIN (low threshold) would be able to rule out the low-mass window of the (g − 2) µ parameter space.
Our results show very promising prospects for separating a potential (g − 2) µ signal of a U (1) Lµ−Lτ from a U (1) Lµ boson through the combination of future data from NA64µ, CEνNS experiments at spallation sources and, crucially, direct dark matter detectors. This strategy would allow us to confirm U (1) Lµ−Lτ as the solution to the (g − 2) µ excess.

A Discovery Lines for LZ and XENONnT
In this appendix, we show the 5σ discovery lines for LZ and XENONnT in Fig. 13 and Fig. 14, respectively. Please refer to Sec. 5 for details. Fig. 13. The same as in Fig. 11 Fig. 11 but for XENONnT.