Solar neutrino probes of the muon anomalous magnetic moment in the gauged $U(1)_{L_\mu-L_\tau}$

Models of gauged $U(1)_{L_\mu-L_\tau}$ can provide a solution to the long-standing discrepancy between the theoretical prediction for the muon anomalous magnetic moment and its measured value. The extra contribution is due to a new light vector mediator, which also helps to alleviate an existing tension in the determination of the Hubble parameter. In this article, we explore ways to probe this solution via the scattering of solar neutrinos with electrons and nuclei in a range of experiments and considering high and low solar metallicity scenarios. In particular, we reevaluate Borexino constraints on neutrino-electron scattering, finding them to be more stringent than previously reported, and already excluding a part of the $(g-2)_\mu$ explanation with mediator masses smaller than $2\times10^{-2}$ GeV. We then show that future direct dark matter detectors will be able to probe most of the remaining solution. Due to its large exposure, LZ will explore regions with mediator masses up to $5\times10^{-2}$ GeV and DARWIN will be able to extend the search beyond $10^{-1}$ GeV, thereby covering most of the area compatible with $(g-2)_\mu$. For completeness, we have also computed the constraints derived from the recent XENON1T electron recoil search and from the CENNS-10 LAr detector, showing that none of them excludes new areas of the parameter space. Should the excess in the muon anomalous magnetic moment be confirmed, our work suggests that direct detection experiments could provide crucial information with which to test the $U(1)_{L_\mu-L_\tau}$ solution, complementary to efforts in neutrino experiments and accelerators.

Models of gauged U (1)L µ−Lτ can provide a solution to the long-standing discrepancy between the theoretical prediction for the muon anomalous magnetic moment and its measured value. The extra contribution is due to a new light vector mediator, which also helps to alleviate an existing tension in the determination of the Hubble parameter. In this article, we explore ways to probe this solution via the scattering of solar neutrinos with electrons and nuclei in a range of experiments and considering high and low solar metallicity scenarios. In particular, we reevaluate Borexino constraints on neutrino-electron scattering, finding them to be more stringent than previously reported, and already excluding a part of the (g − 2)µ explanation with mediator masses smaller than 2 × 10 −2 GeV. We then show that future direct dark matter detectors will be able to probe most of the remaining solution. Due to its large exposure, LZ will explore regions with mediator masses up to 5×10 −2 GeV and DARWIN will be able to extend the search beyond 10 −1 GeV, thereby covering most of the area compatible with (g − 2)µ. For completeness, we have also computed the constraints derived from the recent XENON1T electron recoil search and from the CENNS-10 LAr detector, showing that none of them excludes new areas of the parameter space. Should the excess in the muon anomalous magnetic moment be confirmed, our work suggests that direct detection experiments could provide crucial information with which to test the U (1)L µ −Lτ solution, complementary to efforts in neutrino experiments and accelerators.

I. INTRODUCTION
With the completion of the Standard Model (SM) after the Higgs boson discovery and the absence of any hint for new particles at the energy frontier, the focus of particle physics is shifting towards the exploration of the intensity frontier. At present, increasing emphasis is being put on the emerging discrepancies in precision studies of low-energy observables. In this context, the anomalous magnetic moment of the muon, (g − 2) µ , stands out as one of the most precisely measured quantities in particle physics. Its determination with the E821 experiment at BNL [1][2][3] created a now almost two-decade-old puzzle, as the experimental value deviates from its SM prediction by up to ∼ 3.7 σ 1 [6][7][8]. This has been interpreted as a possible signature of physics beyond the Standard Model (BSM), leading to a plethora of constructions with new physics in the leptonic sector [9][10][11][12][13].
In an effort to shed new light on the true nature of the observed discrepancy, the upcoming measurement at the E989 experiment at Fermilab is projected to determine (g − 2) µ with a precision of 140 parts per billion [14]. This would result in a fourfold improvement over the present experimental error [15]. Furthermore, an independent measurement of this observable at J-PARC [16] with comparable precision to the previous determination at BNL will be able to provide further insight on this potential discrepancy. With these experimental efforts underway and the first data release from E989 imminent, it seems timely to explore complementary probes of the conjectured explanation of the observed shift in (g − 2) µ .
In this article, we focus on a particularly interesting solution proposed in terms of a new light mediator associated with a new U (1) Lµ−Lτ gauge symmetry, which received a lot of attention over the last years [17][18][19][20][21][22][23]. Since gauging the difference between two lepton-flavour numbers is anomaly-free within the SM, these constructions provide a minimal extension of the SM gauge group without the need for extra fermions [24,25]. It has been noted that such a model can simultaneously account for the observed discrepancy in (g −2) µ and even accommodate dark matter with the correct relic abundance [18,[26][27][28][29][30][31]. Furthermore, as suggested in Ref. [32], a U (1) Lµ−Lτ gauge boson A can significantly alleviate the recent 3 σ deviation of local measurements [33][34][35] of the Hubble parameter, H 0 , from the value inferred from early-time cosmology in CMB data [36]. 2 In particular, a modification of early-time cosmology with the presence of an effective component of dark radiation of ∆N eff ∼ 0.4 neutrino-equivalent species could reconcile the two measurements [38].
The new U (1) Lµ−Lτ gauge boson does not have couplings to electrons or quarks at leading order. This makes it hard to test the model in collider or fixed target experiments, which employ electron or hadron beams. However, its tree-level couplings to mu-and tau-flavoured neutrinos open up a unique opportunity to search for such a boson in processes probing neutrino interactions. Almost all of the dominant current bounds on MeV-mass U (1) Lµ−Lτ bosons result from neutrino probes like neutrino trident production [39], neutrino-electron scattering in Borexino [40], neutrino cooling of white dwarfs [19] or N eff during big bang nucleosynthesis [32,41,42]. The only competitive constraints not relying on neutrinos stem from four-muon searches at BaBar [43] and CMS [44].
The aim of our work is to explore the sensitivity of solar neutrino scattering off electrons and nuclei as a future independent probe of the region of parameter space relevant for (g − 2) µ . We start our analysis by considering previous limits derived from Borexino [45,46] data, and show that these have been generally underestimated. We update these bounds, explicitly taking into account the ambiguity arising around the so-called solar metallicity problem, and we perform all computations with solar fluxes as predicted in the solar standard model (SSM) for both the scenario of a low and high metallicity Sun [47].
Next, we consider coherent elastic neutrino-nucleus scattering (CEνNS). This elusive process was first detected at the COHERENT experiment in CsI[Na] nuclei [48,49], an observation which has been recently extended to argon nuclei [50]. These results can be used to constrain the parameter space of new light mediators in the neutrino sector and, in fact, the data from the CsI run was used to derive limits for the U (1) Lµ−Lτ model [51]. In our work we extend the analysis to derive new bounds from the recent LAr results.
Finally, we investigate the potential of direct dark matter detection experiments to probe the region of the U (1) Lµ−Lτ parameter space consistent with the muon anomalous magnetic moment. Although originally designed to search for dark matter particles, direct detection experiments can be sensitive to both neutrino-electron elastic scattering and coherent neutrino-nucleus scattering (the latter is usually considered to be an irreducible background for dark matter searches). In the energy range in which we are interested, the main contribution comes from solar neutrinos. In our analysis we derive constraints from the recent results on electron recoils in XENON1T [52] and argue that they do not explore new areas of the parameter space of the U (1) Lµ−Lτ model. We then investigate the potential reach of upcoming experiments based on germanium and xenon, inspired by SuperCDMS [53] and LUX-Zeplin (LZ) [54], as well as third generation detectors inspired on liquid noble gases, such as DarkSide-20k [55] and DARWIN [56]. We find that the optimal strategy to probe the (g − 2) µ solution is a large exposure, which favours liquid xenon detectors.
This article is organised as follows. In Section II we introduce the minimal gauged U (1) Lµ−Lτ model, we discuss its impact on the muon anomalous magnetic moment, and we compute the contribution to non-standard neutrino interactions. In Section III we describe the various ways to probe this model via neutrino interactions, in particular via neutrino-electron scattering in Borexino (Section III.1), via coherent neutrino-nucleus scattering in the COHERENT experiment (Section III.2), and via both of these processes in direct dark matter detectors (Section III.3). In Section IV we give a detailed discussion of our results, before we conclude in Section V.

II. MINIMAL GAUGED U (1)L µ−Lτ
A peculiar feature of the SM Lagrangian is its invariance under the accidental global symmetries U (1) B,Le,Lµ,Lτ . These global symmetries, if combined into the groups U (1) B−L and U (1) Li−Lj (with i, j = e, µ, τ ) or linear combinations thereof, can be promoted to anomaly-free gauge symmetries without the addition of any extra new fermion fields 3 . Of these extra anomaly-free gauge symmetries the case of U (1) Lµ−Lτ is particularly interesting as it still allows for an explanation of the (g − 2) µ anomaly [17,18].
We begin our discussion of models of extra gauged U (1) Lµ−Lτ symmetries by introducing their generic form of the Lagrangian in the gauge basis, where X α denotes the U (1) Lµ−Lτ gauge boson and B αβ and X αβ are the field strength tensors of the hypercharge U (1) Y and the new U (1) Lµ−Lτ , respectively. The parameters Y and g µτ denote the kinetic mixing parameter and gauge coupling constant. The gauge current of the new symmetry in the minimal setup is given by where L 2 = (ν µL , µ L ) T and L 3 = (ν τ L , τ L ) T denote the second and third generation SU (2) L lepton doublets.
In order to study the phenomenology of a gauged U (1) Lµ−Lτ symmetry, we need to determine the interactions of the mass eigenstate A α of the associated gauge boson X α . The kinetic terms in Eq. (1) are not diagonal and must be diagonalised by a field redefinition of B α and X α . Subsequent rotation to the mass basis 4 yields interaction terms of the following form [19], with the coupling matrix expanded to leading order in the small kinetic mixing parameter and mass ratio of the neutral bosons, Here, θ W denotes the weak mixing angle, M X is the bare mass of the new gauge boson and M SM Z = g 2 + g 2 v/2 is the SM mass of the Z-boson. From Eqs. (3) and (4) we see that interactions of the massless photon state A remain the same as in the SM. The massive hidden photon A , however, couples both to the L µ − L τ current J µ−τ α via gauge interactions and to the electromagnetic current J EM α via kinetic mixing. Hence, the interactions of the physical hidden photon A can generically be expressed as where the coupling coefficients c f for the different fermions f are summarised in Table I. In a simple U (1) Lµ−Lτ extension of the SM, kinetic mixing can arise at tree level through the gauge invariant, renormalisable mixing term in Eq. (1). However, if the new gauge symmetry U (1) Lµ−Lτ is embedded into a larger non-abelian symmetry structure in the UV (such as e.g. SU (N ) Lµ−Lτ → U (1) Lµ−Lτ ) a fundamental mixing term is forbidden by gauge invariance. Nevertheless, kinetic mixing will still be generated radiatively at the loop level and is hence unavoidable. For a general in-detail analysis of the origin of one-loop kinetic mixing we refer to the discussion in Appendix A. For the sake of clarity, we will restrict ourselves in the present work to the case where a tree-level mixing term is either forbidden by the larger symmetry structure of the UV theory, or is subdominant compared to the loop-induced mixing.
At low energies, q 2 M 2 Z , the kinetic mixing can be approximately considered to be only with the photon. In this regime, the full one-loop result for the loop-induced kinetic mixing parameter in the minimal U (1) Lµ−Lτ setup reads, where the fermions µ and τ , charged under both U (1) Lµ−Lτ and the SM U (1) EM , are running in the loop. Here q 2 denotes the momentum transfer flowing through the loop. For very large momenta, q 2 m 2 τ , the mixing is momentum suppressed, µτ ∝ m 2 µ /q 2 − m 2 τ /q 2 . However, as long as the momentum transfer is below the muon mass q 2 m 2 µ , the mixing is approximately constant 5 , which corresponds to the physical situation encountered in all the processes discussed in this work. 5 It should be noted that the kinetic mixing parameter as defined in Eq. (8) is negative, µτ < 0.

II.1. Muon anomalous magnetic moment
Due to its gauge interactions with the second generation leptons the U (1) Lµ−Lτ gauge boson A contributes to the anomalous magnetic moment of the muon, a µ = (g − 2) µ /2, via the loop process displayed in Fig. 1. For any neutral gauge boson with vectorial couplings to muons (as in U (1) Lµ−Lτ ), the additional contribution to a µ can be expressed in the compact form [58,59], where α = g 2 µτ /4π, x µ = m µ /M A and Q µ denotes the charge of the muon under the new gauge symmetry, i.e. in the case of U (1) Lµ−Lτ we have Q µ = 1.
These values expose the quoted 3.7 σ deviation [8] of Interpreting this deviation as a signal for new physics, the observed discrepancy can be explained by virtue of Eq. (9) with the presence of a U (1) Lµ−Lτ gauge boson for an appropriate choice of the gauge coupling strength g µτ . It is worth noting that the additional contribution ∆a µ induced by vectorial couplings in Eq. (9) is strictly positive, regardless of the charge of the muon Q µ under the new symmetry. Hence, models with vectorial couplings can only explain an upward fluctuation in (g − 2) µ . With new experimental results from E989 imminent it is worthwhile pointing out that in case of a possible observation of a downward fluctuation compared to the SM prediction, models with vectorial couplings to muons, like U (1) Lµ−Lτ , would be in strong tension with the experimental findings. In this sense, (g − 2) µ could potentially be used to rule out such models.
In case the central value of the observed anomalous magnetic moment should move downward and would be in agreement with the SM prediction, (g − 2) µ would simply turn into a constraint and yield an upper bound on the gauge coupling constant g µτ of the U (1) Lµ−Lτ boson.

II.2. Non-standard neutrino interactions
The leptophilic U (1) Lµ−Lτ gauge boson induces new neutrino-matter interactions. If the mass of the new mediator A is much heavier than the relevant energy scale of the scattering, M 2 A q 2 , 6 For details on the individual contributions entering the SM calculations we refer the reader to Ref. [8].
the mediator can be integrated out and the process be described by a four-fermion contact interaction. However, even in the regime where the mediator mass and the relevant scattering energy are comparable, M 2 A q 2 , we can encapsulate the hidden photon interaction in a quasi-effective (energy-dependent) contact term of the neutrinos ν α and fermions f . Making use of the framework of neutrino non-standard interactions (NSI) [78,79] this contact term can be written as 7 , where G F denotes the Fermi constant and P = {P L , P R }. In order to match the gauge boson interaction onto that of the effective operator in Eq. (13), we need to identify the amplitudes of the [ν α γ ρ P L ν α ] f γ ρ P f interaction in the effective theory with the one in the full theory. Making use of the expression for the four-momentum transfer, q 2 = −2 m f E R , the matching condition yields an effective energy-dependent NSI coupling of , for a kinetic mixing coupling to f , where the minus sign in the kinetic mixing coupling is due to the relative sign difference of the gauge and kinetic mixing interaction of the A in Eqs. (3) and (4). Furthermore, we denote the corresponding neutrino charges under U (1) Lµ−Lτ by Q να = {0, 1, −1} for α = {e, µ, τ }, respectively. In the following we will use the NSI framework to study the effect of an additional U (1) Lµ−Lτ gauge boson on a number of neutrino physics processes.
• Coherent elastic neutrino-nucleus scattering The first neutrino interaction we are considering in this article that receives observable modifications in the presence of a U (1) Lµ−Lτ gauge boson is CEνNS. In order to compute the A contribution to this process we can make use of the NSI formulation in Eq. (13). This modifies the low-energy formulation of the weak neutral currents to where we denote the coupling of the SM Z-boson to the fermion f as Following the standard derivation of the cross section [80,81], we can write this parametrically as where E ν is the energy of the incoming neutrino, E R and M N are the recoil energy and mass of the scattered nucleus. For simplicity, we have assumed that the electromagnetic form factors F Z (q 2 ) and F N (q 2 ) describing the charge distribution within protons and neutrons are roughly given by the Helm nuclear form factor F (E R ) [82,83], which describes the nucleon distribution within a nucleus in terms of the magnitude of the three-momentum transfer 7 Note that here ε f P αα denotes the strength of the effective NSI, not to be confused with the loop-induced kinetic mixing parameter µτ (q 2 ).
Tree-level contributions to the scattering ναe → ναe. The first diagram only contributes to the process for α = e, whereas the last one contributes only for the flavours α = µ, τ .
where j 1 denotes the spherical Bessel function of order one, and R 2 0 = c 2 + 7 3 π 2 a 2 − 5s 2 with c (1.23 A 1/3 − 0.6) fm, a 0.52 fm and s 0.9 fm. 8 The relevant neutrino-nucleus scattering amplitudes corresponding to SM and NSI scattering, respectively, can be expressed as where we have made use of the vector couplings defined as Replacing our expressions for the A -induced NSI couplings defined in Eq. (14), the full expression for the differential neutrino-nucleus scattering cross section can be written in the compact form, where we have defined the coherence factor for the effective neutrino-nucleus coupling via the SM Z-boson as From Eq. (23) we can see that the U (1) Lµ−Lτ gauge boson couples proportional to the total electric charge of the nucleus, Z. This is expected as in the case of U (1) Lµ−Lτ the coupling to hadrons is generated from kinetic mixing and hence proportional to the photon coupling, i.e. the electric charge 9 .
• Elastic neutrino-electron scattering The second process we consider to probe neutrino interactions of a new U (1) Lµ−Lτ is elastic neutrino-electron scattering, ν α e → ν α e. In the SM, neutrino-electron scattering proceeds via the weak interactions. For the neutrino flavours α = µ, τ this process is solely mediated by the neutral current interaction in Eq. (15) (corresponding to the diagrams in the centre and right panel of Fig. 2). However, for α = e this process is also mediated by a charged current interaction (corresponding to the diagram in the left panel of Fig. 2). The corresponding low-energy charged current contact term reads The general form of the neutrino-electron scattering cross section induced by these low-energy neutral and charged current interactions can then be expressed as [85][86][87][88][89][90], where once more E ν denotes the initial energy of the incoming neutrino and E R the recoil energy of the scattered electron. In the SM the relevant weak interaction couplings are given by for α = µ, τ , where g e L and g e R refer to the couplings of the Z-boson to electrons as defined in Eq. (16). It is now straightforward to extend Eq. (26) to incorporate neutrino-electron scattering mediated by the new U (1) Lµ−Lτ . By virtue of the NSI formulation of the A interaction, we can simply make the replacement where ε eP αα are the A -induced NSI couplings to electrons defined in Eq. (14). As the first generation of leptons is uncharged under U (1) Lµ−Lτ the cross section for ν e e → ν e e scattering will to lowest order remain unchanged. However, for α = µ, τ the scattering cross section for the process ν α e → ν α e is modified and explicitly reads The first term in the curly braces 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 να of the scattered neutrino ν α .
Apart from the two types of elastic neutrino scattering processes just discussed, quite universal limits on NSI couplings are typically derived from the induced neutrino-matter potentials, which can be probed e.g. in neutrino oscillation experiments [91,92]. These constraints are for example quite stringent on hidden photon models featuring gauge couplings to baryons [93][94][95]. However, in the case of U (1) Lµ−Lτ the hidden photon couplings to matter (i.e. p, n and e) arise only via kinetic mixing and are therefore proportional to the electric charge. Hence, as the matter within the Earth or experimental apparatuses is electrically neutral, the net effect of NSI couplings proportional to the electric charge is exactly zero as they cancel out [96], Therefore, the typically strong limits from NSI matter potentials do not apply in the case of a U (1) Lµ−Lτ gauge boson. Nevertheless, NSIs generated in the U (1) Lµ−Lτ model still have the potential to affect global fits of neutrino mixing parameters, as they directly change the neutrino-electron elastic scattering cross section. This is important, as some of the main results used to constrain the mixing parameters θ 12 and ∆m 12 are the measurements of the scattering of 8 B neutrinos with electrons by the Super-Kamiokande, SNO+ and Borexino experiments [97][98][99]. Any new physics which affects the electronneutrino scattering rate will also affect this result [100], leading to an incorrect determination of these neutrino mixing parameters. This kind of degeneracy has been explored in the context of quark NSIs [92], but it has the potential to be much more significant with electron NSIs.
Two potential ways to resolve this degeneracy are apparent. Firstly, the effects of changing the neutrino mixing parameters and of introducing non-zero NSIs are likely to be different for measurements taken in different energy windows. All of the 8 B neutrinos measured at Super-Kamiokande had energies above 5 MeV. This means that they had passed through the MSW resonance and so the probabilities of finding them in each neutrino flavour eigenstate will be different from those for lower energy solar neutrinos, while an energy dependence is also found in the neutrino-electron cross section. Secondly, for the case of a U (1) Lµ−Lτ model the NSIs generated with electron-neutrinos are significantly suppressed relative to muon and tau neutrinos. The other major result used in calculating the quoted neutrino mixing parameters is the one from KamLAND, which detects reactor neutrinos by inverse-beta decay [101]. This means it is only able to detect electron anti-neutrinos, and is therefore unaffected by the NSIs e µµ and e τ τ . In this context, it is most interesting to note that a tension currently exists between the mixing parameters computed using only KamLAND data and those using only solar neutrino data, especially in the value of ∆m 12 [102]. It seems likely that the only consistent way to obtain a simultaneously correct treatment of both neutrino mixing and NSIs with quarks and electrons, suitable to determine whether NSIs can resolve the tension in measurements of ∆m 12 , would be to include all the mixing parameters and NSI couplings in a single global fit analysis. Once such a holistic global fit including all NSI parameters will be available, it will allow for an improved analysis of the limits. In the absence of such a fit and for the remainder of this work we take our neutrino mixing parameters from the global fit in Ref. [102].

III. SOLAR NEUTRINO PROBES OF THE GAUGED U (1)L µ−Lτ
As we have seen in the previous section, the U (1) Lµ−Lτ hidden photon introduces new neutrino interactions that have an effect on neutrino electron-scattering and coherent neutrino-nucleus scattering. In this section we will study how these effects can be probed. We start by considering neutrino-electron scattering in dedicated neutrino experiments and reevaluate constraints from the Borexino detector. Then, we consider coherent neutrino-nucleus scattering in spallation facilities and derive new limits from the recent CENNS-10 LAr results. Finally, we address dark matter experiments sensitive to both electron and coherent-nucleus scattering. In this context, we first derive limits from the recent XENON1T electron recoil data before we study the prospects for a number of upcoming and future direct detection experiments.
These three types of experimental techniques cover the parameter space of the gauged U (1) Lµ−Lτ model in complementary ways. We will now explore how effectively they can probe the area of the parameter space that is consistent with a solution to the muon anomalous magnetic moment.
In the case of Borexino and the various dark matter experiments we will consider, the most relevant source of neutrinos is the Sun. Solar neutrinos can be divided into several individual flux spectra, each associated with a different stage in the fusion process. The most important fluxes for us are those generated during the proton-proton (pp) chain: the pp, 7 Be, pep, 8 B and hep neutrinos. Each of these spectra has a different total flux, which must be properly modelled for us to be able to place reliable constraints. The values we use for the individual neutrino fluxes are computed in Ref. [47] within the context of the Standard Solar Model (SSM).
A significant source of uncertainty in the SSM prediction is the as-yet unsolved solar metallicity problem [103][104][105][106][107]. The rates of various fusion processes depend on the metal content 10 in the Sun. Two possible scenarios exist: the so-called high metallicity (HZ) GS98, and low metallicity (LZ) AGSS09 models [47]. The solar metallicity has a significant effect on the rate of the Carbon-Nitrogen-Oxygen (CNO) fusion cycle: the associated CNO neutrino fluxes are around 50% larger in the HZ scenario [47], and a direct measurement of these as-yet unobserved neutrinos could help to finally resolve the solar metallicity problem. However, it also affects the various pp-chain fluxes by up to 10%. The reach of future experiments to explore new physics will therefore differ in these two scenarios. Constraints obtained from current results such as Borexino can differ even more: what appears like a small upwards fluctuation above the SM result under the assumption of a HZ Sun can be a much larger enhancement if the Sun has low metallicity, which could allow for a larger contribution to the interaction rate from new physics.
In order to maintain full generality when deriving limits from solar neutrino scattering we will hence explicitly perform our calculations for both the scenario of a high and low meallicity Sun.

III.1. Borexino
With its large fiducial volume and low background rates, the Borexino experiment has been responsible for many of the most precise direct measurements to-date of the various solar neutrino fluxes through their scattering with electrons in a scintillating liquid target. With analyses covering energies from 0.19 to 20 MeV, it has measured or placed constraints on the fluxes of pp, 7 Be, pep, 8 B, hep, and CNO neutrinos [97,108], leading to constraints on models of new physics which affect the neutrino-electron scattering rate.
Among the various sources of solar neutrinos, Borexino has measured the 7 Be rate to the highest level of precision. This, combined with the fact that the greatest change to the cross section from a new light mediator is at low energy, makes 7 Be neutrinos the best candidate for placing constraints on the gauged U (1) Lµ−Lτ model 11 .
In Ref. [109], the Borexino collaboration measured the 7 Be flux with a precision of 5%, finding it to be in agreement with the SM and SSM. This result has been used in an earlier work to place constraints on a model of gauged U (1) B−L , excluding all regions of parameter space where the enhancement of the flux over the SM prediction exceeds 8% (the corresponding 90% CL given a 1σ uncertainty of 5%) [45]. This result has itself been used to place a constraint on a U (1) Lµ−Lτ model by remapping the associated couplings [40]. We have improved on this method in several ways. Firstly, by calculating constraints directly for a U (1) Lµ−Lτ model we are able to take into account the full scattering cross section from Eq. (29) including the interference term. Secondly, the original B − L result did not take into account the uncertainty on the theoretical prediction from the SM and SSM, including the uncertainty stemming from the solar metallicity problem. Finally, results from Phase II of Borexino have since measured the flux with even higher precision, of around 2.7% [108]. We calculate constraints on a U (1) Lµ−Lτ model from both Phase I and II Borexino results in both metallicity scenarios using a chi-squared test. In Appendix C the same process is used to derive updated constraints on the U (1) B−L model for which these constraints have been first derived.
Under the assumption of high solar metallicity, and based on the work of [47], the Borexino collaboration computed a theoretical uncertainty on the 7 Be rate of 5.8%. Combining this with the uncertainty on the experimental measurements, we derive constraints at the 90% confidence level using a chi squared test. While the two Borexino results are consistent with each other and with the SM (under the assumption of a high metallicity SSM), the Borexino Phase II measured a higher value for the 7 Be flux, and so gives weaker constraints on our model. As these two results can be consistently explained as a downwards and upwards fluctuation, respectively, we will display the constraints derived from each analysis, with the assumption that a more sophisticated analysis would lie somewhere between the two. The predicted 7 Be neutrino flux is ∼ 10% lower under the assumption of a low metallicity Sun than it is with a high metallicity Sun. A larger contribution to neutrino-electron scattering from new physics is therefore allowed by the measurements from Borexino.
A more complete method of deriving constraints on a U (1) Lµ−Lτ model from all the solar neutrino measurements from Borexino would be to combine them all in a single chi-squared analysis. However, as the best fit values for the count rates of pp, 7 Be, and pep were computed simultaneously in a single multivariate fit [108], a careful handling of the various uncertainties would be required, and we leave this to future works.

III.2. Coherent neutrino-nucleus scattering at COHERENT
Having considered dedicated searches for neutrino-electron scattering, we turn our attention to experiments which are specifically designed to measure coherent neutrino-nucleus scattering. The CO-HERENT collaboration observed this elusive process for the first time in CsI [49] and, more recently, measured it again using a liquid argon target (in the CENNS-10 LAr detector [50]) at the Oak Ridge National Laboratory Spallation Neutron Source. The results are compatible with the SM prediction and therefore can be interpreted as constraints on any new physics contribution, which are especially relevant for models with light mediators [51,[110][111][112][113]. For completeness, we have incorporated the limits from the CsI run, derived in Ref. [51]. To derive the constraints from the new LAr data, we proceed as follows.
The number of expected CEνNS events can be expressed as where dN να /dE ν refers to the spectrum of the incoming neutrino ν α produced at SNS, the minimum and maximum neutrino energies are E min ν ≈ M N E R /2 and E max ν = m µ /2, and the expression for the differential cross-section is that of Eq. (23). The number of target nuclei is N targ = M tot /M N , with M tot = 24 kg (the exposure of CENNS-10 LAr) and M N the atomic mass of 40 Ar, which we consider (for simplicity) to have 100% isotopic abundance. For the energy-dependent efficiency, (E R ), we take the one given in analysis A of Ref. [50].
The spectrum of stopped-pion neutrinos produced at SNS contains a delayed neutrino flux, composed of muon antineutrinos and electron neutrinos, both with a maximum neutrino energy of E max ν = m µ /2 (see e.g., Refs. [49,114]) and fluxes given by as well as a monoenergetic muon neutrino flux, with E ν = 29.8 MeV, In these expressions, 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) (the 6.12 GWh exposure corresponds to N POT = 1.37 × 10 23 ). The number of neutrinos produced per POT is given by r = 0.08, while L = 27.5 m is the CENNS-10 baseline. In order to interpret the experimental data we need to translate from nuclear recoil energies to the equivalent energies of electron recoils. The quenching factor allowing for this conversion from nuclear recoil energies to electron-equivalent energies is parametrised as Q F = 0.246 + 7.8 × 10 −4 E R , such that E[keV ee ] = Q F E R [50].
In order to derive the corresponding limit from the CENNS-10 LAr run on the parameter space of a U (1) Lµ−Lτ gauge boson we perform a chi squared test. For each point in the (M A , g µτ ) parameter space, the χ 2 is then computed as follows [113,115] where the number of observed events in the detector is N meas = 159 ± 43 (stat.)±14 (syst.), the number of background events is N bgd = 563, and σ α = 0.085 accounts for a scaling systematic.

III.3. Direct dark matter detection experiments
Dark matter direct detection experiments are, in principle, also sensitive to MeV-scale neutrino interactions through both their scattering off electrons and their coherent scattering off nuclei [116][117][118]. Although their performance is still limited by the somewhat smaller target size than dedicated neutrino experiments, the next generation of detectors with extremely low energy threshold (such as SuperCDMS) or increased fiducial volume (liquid noble gas detectors such as LZ or DarkSide) will provide a complementary test of new physics in the neutrino sector.
The original focus of most of the detectors was to search for nuclear recoils (NR) for which the experimental design and data analysis generally guarantees a very small (ideally zero) background thanks to a very careful rejection of the large electron recoil background. However, in the last decade a great effort has been made to exploit electron recoils (ER) for dark matter (and other new physics) searches. Electron recoils allow to probe lighter particles but, in general, at the expense of introducing a larger background. In our analysis we consider both types of searches, as they explore the parameter space of new physics in the neutrino sector in a synergic manner [116]. In fact, when we give up on discrimination between both types of events, one can do a combined (NR+ER) analysis in which the signals (and backgrounds) from both types of recoils are added. This can be advantageous to improve the threshold for nuclear recoils, so we also take it into account.
In full analogy to Eq. (31) we write the expected number of solar neutrino-electron or -nucleus scattering events at direct detection experiments as where now E min ν = (E R + E 2 R + 2m T E R )/2 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, n T is the total number of targets (electrons or nuclei) per unit mass 12 , and dσ να T /dE R is the differential cross section of the elastic scattering of neutrinos off electrons (Eq. (29)) or nuclei (Eq. (23)). We sum across all neutrino flavours, α = e, µ, τ , and the incoming (solar) neutrino flux, dφ νe /dE ν , is multiplied by the transition probability between the source and the detector P (ν e → ν α ). For simplicity, we will consider the experimental exposure, ε, to be energy-independent. The threshold energy, E th , and maximum energy E max are specific to the experiment.
So far, none of the results published by direct detection experiments probe new areas of the parameter space of the U (1) Lµ−Lτ model. For this reason, in our analysis we have decided to focus on the potential of some representative experiments which are currently being deployed (SuperCDMS and LZ) and on third generation proposals (DarkSide-20k and DARWIN). However, the most sensitive search for electron recoils has been recently released by the XENON1T collaboration [52], which seems to hint at a possible excess at low energies. Due to its relevance and timeliness, we have decided to include this data in our analysis as well.

III.3.1. XENON1T
The XENON1T collaboration [119] operated a liquid xenon time projection chamber (with an active target of 2 tonnes) from 2016 to 2018 to look for exotic events at the Laboratori Nazionali del Gran Sasso (LNGS). Particle interactions within the detector produce a prompt scintillation (S1) signal, accompanied by a delayed electroluminescence (S2) signal, both detected by an array of photomultipliers. The ratio of both signals (S2/S1) can be used to discriminate electron from nuclear recoils, thereby guaranteeing a low background for the latter in the canonical analysis mode (which regards electron recoils as background).
In a recent analysis [52] the data for electron recoils in XENON1T Science Run 1 has been analysed to look for potential contributions from new physics. The greatest challenge of this study is the existence of an irreducible background dominated at low energies by the β decays of 214 Pb and with additional contributions from other sources of contamination as well as the xenon target itself. The different contributions to the background are determined using a global fit in the region 1 − 210 keV to obtain a best fit model, denoted as B 0 .
The results obtained by the XENON1T collaboration are shown in Fig. 3, together with the best fit model for the background (in red). An excess is observed in the low energy region (below 7 keV and more prominent in the first two bins) which can be understood as a hint for new physics (for example, solar axions, a neutrino magnetic moment, or light bosonic dark matter). However, as the authors of [52] point out, a plausible explanation is also a previously unconsidered source of background, caused by traces of tritium (which undergoes β decay) in the detector.
The enhancement of the spectrum over background attributed to an additional U (1) Lµ−Lτ boson in our region of interest of the parameter space would be mostly flat over the energy range explored by XENON1T. In Fig. 3 we choose two benchmark points to show the enhancement to the recoil spectrum produced by a U (1) Lµ−Lτ model over the top of the best-fit background spectrum from Ref. [52], B 0 . For the first benchmark point, BP1, we take M A = 15 MeV and g µτ = 5.6 × 10 −4 , which lies in the region which can simultaneously resolve both the H 0 and (g − 2) µ anomalies, and is not excluded by any previous bounds. For the second, BP2, we take the same mass but increase the coupling to g µτ = 1.7 × 10 −3 , which is the point at which our model disagrees with the data at the 90% confidence level (according to the calculation described below). This means that a U (1) Lµ−Lτ boson cannot explain the observed fluctuation in the XENON1T data, which is peaking in the region between 2 − 3 keV. A U (1) Lµ−Lτ or other light vector mediator could produce a feature resembling the XENON1T data, but this would require a much smaller mass, below 100 keV [116], a region of the parameter space which is already well-constrained in our model [32,120,121].
Therefore, we use the data from XENON1T to derive an approximate upper limit on the parameter space of the U (1) Lµ−Lτ boson. Since our spectrum is relatively flat and featureless, we perform a simple unbinned delta chi-squared (∆χ 2 ) analysis. Summing the data points, we compute the total number of events observed in the energy range 2 − 30 keV, N tot . We then compare this value with the number of events predicted in the B 0 background model, and in the B 0 +U (1) Lµ−Lτ scenario. We then place an upper limit at the 90% confidence level using a ∆χ 2 test.
For the purpose of performing the chi squared test we use the number of background events as determined in the background-only fit by the Xenon collaboration. However, it should be noted that there is a large systematic uncertainty on the number of expected background events coming from β decays of 214 Pb, quoted as N exp Pb = (3450, 8530). This is particularly interesting as this constitutes by far the most dominant background in the low energy part of the spectrum, where e.g. the number of background events from SM solar neutrinos N exp Sν = 220.7 ± 6.6 is more than an order of magnitude smaller. Due to the poor knowledge of the lead background it is fitted to data in the background-only model. Taking this as a fixed background for the purpose of our chi squared test with additional solar neutrino scattering due to a U (1) Lµ−Lτ boson is overly simplistic, as any additional events from hidden photon-induced solar neutrino scattering could be easily compensated for by a fluctuation in the lead background. However, the approximation we take will lead to an overestimate of the severity of the constraint, and since the limit we derive does not constrain any new regions of the parameter space we can be confident that any limits derived from this XENON1T result will not be competitive with existing constraints. We also perform an equivalent analysis of a U (1) B−L model, and as in the U (1) Lµ−Lτ case we find that limits from this result are not competitive (see Fig. 8).
In a more careful derivation of the limit from the XENON1T result one should perform a full profile likelihood analysis, allowing the backgrounds to vary within the expected range as nuisance parameters. However, such a dedicated statistical analysis of the limits is beyond the scope of this sensitivity study.

III.3.2. Forthcoming and future detectors
Let us now turn our attention to the future sensitivities of dark matter experiments. We consider simplified setups based on the characteristics of some prominent current detectors (to which we will refer as second generation, or G2) and extend the analysis to some next generation (G3) proposals. The experimental configurations are shown in Table II, where we specify the detector target, exposure, and energy ranges used in each analysis. Due to the different energy-loss mechanisms associated with each type of recoil event, the signals produced by electron and nuclear recoils must be interpreted differently when translating them into recoil energies. We differentiate between nuclear-recoil energies (denoted by keV nr ) and electron-equivalent energies (denoted by keV ee ).
In general, far more background events are expected from electron recoils (ER) than nuclear recoils (NR). A common way of reducing backgrounds in searches for nuclear recoils is through the combination of multiple detection channels, with nuclear and electron recoils having different profiles when viewed through both simultaneously. However, this typically limits the minimum energy threshold of the analysis. We consider three classes of analyses, denoted in Table II as NR, ER, and NR + ER. In the first two, our energy ranges are chosen to allow maximum discrimination between NR and ER events, and the signals of each recoil type are studied separately. In the latter, the two signals are combined, sacrificing a low-background nuclear recoil analysis to allow us to probe lower energy recoils. In this last type of analysis, we must first convert all recoil energies to the nuclear-recoil equivalent before combining our signals.
The nuclear recoil spectra are integrated up to the maximum kinematically allowed energies for the NR analyses, which follow from a collision with the highest energy solar neutrinos available (here the tail of the hep neutrinos). Finally, we have not included energy resolution effects. The G2-Ge setup is based on the planned configuration of SuperCDMS at SNOLAB [53]. The Super-CDMS experiment includes two types of detectors: iZIPs, which allow for an excellent discrimination between electron and nuclear recoils, and HV [122], which operate at high voltage to exploit phonon amplification and achieve a much lower energy threshold at the expense of not discriminating between electron and nuclear events 13 . The exposure corresponds to five years of operation with an 80% live time and the background levels for electron and nuclear recoils in each detector configuration are taken from the predictions in Ref. [122]. We include an overall 75% (85%) signal efficiency for the iZIP (HV) detectors after analysis cuts [53], though we neglect any energy dependence in the efficiency. In our analysis, we compute separate bounds for electron and nuclear recoils using the iZIP configuration and a joint limit (adding electrons and nuclear recoils) for the HV detectors. The nuclear-recoil threshold energies, both assumed to have been derived from a phonon signal, have been converted into electronequivalent thresholds as per Eq. (8) of Ref. [122] using the Lindhard model for the ionisation yield [123]. The maxima of the energy windows resemble the approximate range of the detectors for the ER and NR+ER analyses.
The G2-Xe and G3-Xe experiments have been inspired by the upcoming multi-ton liquid xenon (LXe) experiments LZ and DARWIN, respectively, and also serve as a proxy for competing experiments XenonNT [124] and PandaX [125] using the same target. Both of their target nuclear-recoil thresholds are above the maximum of the solar neutrino spectra; they have been designed as such to minimise this signal, which is a background for WIMP searches. However, the LUX collaboration has been able to set thresholds as low as 1.1 keV nr , with a 3.3 keV nr energy threshold at 50% detector efficiency in both the S1 and S2 channels required for ER/NR discrimination [126]. In light of this and the advent of new analysis techniques allowing for even lower energy thresholds [127], we have set thresholds for both LZ and DARWIN to be similar to this larger result. Above this threshold we can perform a nuclear recoil analysis with 99.5% rejection of ER backgrounds with a 50% acceptance of NR signal events [54,128]. We take our backgrounds from Ref. [54] for LZ and Ref. [128] for DARWIN. In both cases, we extend our analyses beyond the energy range for which backgrounds are specified, taking the conservative assumption of a flat background spectrum within our energy range [129]. The lowest energy measurements in LXe to date have succeeded in detecting an ionisation (S2) signal at energies equivalent to 0.3 keV nr for nuclear recoils [130] and 0.186 keV ee for electron recoils. Abandoning NR/ER discrimination which requires an S1 scintillation signal, we have set lower thresholds at approximately twice this former value for the G3-Xe setup and at 0.7 keV nr for the G2-Xe configuration -the S2-only threshold achieved by XENON100 [131]. We have conservatively assumed that the thresholds for these experiments for the NR+ER analyses are equal in value in both keV nr and keV ee to avoid extrapolating the Lindhard model to energies at which it has not been tested experimentally 14 . The maximum recoil energy has been set at 30 keV ee (∼ 100 keV nr using the Lindhard model) -the energy at which the double-β decay of 136 Xe, a considerable background for all LXe detectors, is expected to dominate over the solar neutrino signal [54].
Finally, the G3-Ar configuration has been based on the future DarkSide-20k detector, for which we also consider a five-year running time. The energy threshold values required for signal discrimination were found to be too high to give competitive constraints from the nuclear recoil spectrum [132], so only the ER-analysis was performed at this threshold. Since the only benefit of taking this higher threshold would be to minimise the background on an NR analysis, we have also taken the much lower threshold required for an ionisation (S2) signal with no ER/NR discrimination of 0.6 keV nr , as achieved by the DarkSide-50 experiment [132]. The maximum of the energy window have been taken from the experimental design report [55] for the ER analysis and placed at 15 keV nr for the NR+ER analysis. The backgrounds we take from Fig. 7 of Ref. [132]. As these do not extend to our maximum energy, we assume that the spectrum is flat above 15 keV nr , which we justify by comparison with Fig. 3 in the same paper.
For each of these setups, we have derived a 90% C.L. exclusion line assuming that no signal for new physics is observed. For a fixed value of the mediator mass, M A , we determine the value of the coupling, g µτ , for which 90% of hypothetical experiments would expect to see an excess over the SM prediction. We have assumed that N follows a Poisson distribution and varied the mean number of events to calculate this. The resulting exclusion lines for each experimental configuration are shown in Fig. 4. We separate the results for electron and nuclear recoils, and confirm (as previously pointed out in the literature [116]) that electron recoils are optimal to test very light mediator masses, whereas nuclear recoils excel for masses above m A > 0.01 GeV. 15 As we will see in the next Section, in the U (1) Lµ−Lτ scenario, experimental constraints rule out the regions with very light mediator masses (in the coupling range suitable to explain (g − 2) µ ). Therefore, the best limits derived from direct detection are mostly based on the NR (or NR+ER) analysis. In this regime, the most advantageous strategy is 14 We have explicitly checked that extrapolating the Lindhard model to low energies would yield at most a 30% improvement to our limits in the very low mass plateau, where M A 10 −4 GeV, and a negligible difference in the region of interest for our work (i.e. for M A 10 −3 GeV). 15 The sensitivity of direct detection experiments to new physics in the neutrino sector can also be interpreted as an increase in the so-called neutrino floor, below which neutrino-nucleus coherent scattering becomes a dominant background for dark matter searches [84,133].
to increase the experimental exposure, rather than further decreasing the threshold. This favours large experiments (such as LZ or DARWIN) versus smaller low-threshold detectors.

IV. DISCUSSION
After having detailed in Section III how neutrino-electron and -nucleus scattering can be used at various different experiments to constrain models of gauged U (1) Lµ−Lτ , we want to discuss our results in the following section.
We show existing bounds on the U (1) Lµ−Lτ gauge boson from white dwarf cooling, neutrino trident production, coherent scattering and four-muon searches in grey in Fig. 5. Concerning neutrino trident production, the most stringent bound is in principle due to the CCFR result [39]. However, it has been pointed out in [134] that the corresponding analysis did not take into account a background from diffractive charm production and hence needs further review. Thus, we show the corresponding limit obtained by CCFR only as a dashed black line in Fig. 5 until this issue is resolved. The region favoured by the observed deviation in (g − 2) µ at the 2 σ level is illustrated by the green band. The region in parameter space explaining the tension in H 0 , corresponding to ∆N eff ∼ 0.2 − 0.5 [32], is shown as the blue band. Furthermore, we take the contour of N eff = 4 derived in [32] as the conservative limit where a model of gauged U (1) Lµ−Lτ should be firmly ruled out by cosmological observations. We display our limit derived from the liquid argon run of COHERENT outlined in Section III.2 as the cyan area. It can be seen that it has a similar sensitivity than the previous caesium-iodide run, but is slightly less constraining. The yellow line represents the limit derived from the recent electron recoil data from XENON1T [52]. As discussed in Section III.3.1 the simplistic treatment of the background in the chi squared test performed to derive the limit leads to an overly optimistic constraint. This limit lies well within previously excluded parameter space, and a more careful treatment of the background would only lead to a weaker bound. We therefore conclude that the current XENON1T constraint is subdominant with respect to the other existing bounds in the U (1) Lµ−Lτ parameter space.
In the top panel of Fig. 5 we show our reevaluated limits from Borexino in red for the case of a high metallicity Sun. The solid red line corresponds to the limit derived with the more recent 2017 Borexino data set [108]. The dashed line corresponds to the updated limit using the older 2011 data set [109], including a more complete treatment of the theoretical and experimental uncertainties. We can see that the limit derived from the more recent data set is less constraining, due to a difference in the two measurements of the scattering rate which is explainable by a statistical fluctuation. Hence, we take the band between the solid and dashed line as the approximate envelope of the true Borexino limit, which should be derived with a dedicated analysis combining the two results in a single chi-squared analysis with appropriate handling of any shared systematic uncertainties. Nevertheless, even in the more conservative case the updated Borexino limit is significantly more constraining than the previously derived bound in [46], which is shown by the red dotted line. In particular, the updated limit rules out previously untested parameter space relevant for the explanation of (g − 2) µ above and hence yields the most stringent low-mass bound on this region of parameter space. Most noticeably, it already excludes part of the region where the (g − 2) µ and H 0 bands overlap and a simultaneous explanation of the two would be in principle possible. The Borexino limits for the case that the Sun is a low metallicity star are shown in the lower panel of Fig. 5. The general picture is quite similar to the high metallicity case. However, as the neutrino fluxes are generally lower in the low metallicity case the corresponding SSM prediction is also lower in the low metallicity case. Overall this leads to less constraining bounds from Borexino, which in the most conservative case rule out U (1) Lµ−Lτ gauge bosons with Nevertheless, Borexino yields constraints which are competitive with those from white dwarf cooling, the most stringent previous bounds on the region of parameter space where (g − 2) µ is still allowed.  Most interestingly, the Borexino bounds still allow for a significantly larger part of the overlap region of the (g − 2) µ and H 0 bands in the low metallicity case. The intriguing fact that a simultaneous explanation of the two tensions is still allowed in both the high and low metallicity scenario makes this region where the two bands overlap a prime target for future experimental searches.
Motivated by the high sensitivity of solar neutrino scattering at Borexino we next want to discuss future sensitivities of various experiments to probe the relevant region in parameter space. In Fig. 6 we show projections for various different experiments on the relevant parameter space. The coloured solid lines show the envelope (i.e. the most constraining bound for a given mass M A ) of the projected limits obtained from analysing electron recoils (ER), nuclear recoils (NR) and the combination of the two (NR+ER) (cf. Fig. 4 for the individual limits) for the various dark matter direct detection experiments discussed in Section III.3.2. The results for SuperCDMS, LZ, DARWIN and DarkSide-20k are illustrated by the red, blue, turquoise and pink lines, respectively. For comparison we show previously derived projections from kaon decays at NA62 [134], neutrino trident production at DUNE [135,136] and a future 10 t · yr run of COHERENT with a NaI/Ar target [51] as black dotted lines.
The projections for the direct detection experiments in the top panel of Fig. 6 were derived using solar neutrino fluxes for the case of a high metallicity Sun. To begin with, it can be seen that SuperCDMS will not be able to improve over the updated Borexino limits. The envelope for SuperCDMS in the shown region is mainly due to the HV analysis not discriminating between electron and nuclear recoils and the iZIP analysis of nuclear recoils only. We explicitly checked and found that the major limitation of SuperCDMS is a nuclear background of 206 Pb decays. Trying to reduce this background further would have the most effect on the limit and potentially allow SuperCDMS to gain sensitivity of unprobed parameter space. Next, we notice that DarkSide-20k is very similar in reach compared to SuperCDMS. The achieved sensitivity is entirely due to the combined NR+ER analysis, as in the setup where discrimination of ER versus NR is possible, the corresponding threshold is too high to see solar neutrino scattering events through nuclear recoils. In contrast to this, we observe that LZ and, most prominently, DARWIN will achieve the best sensitivity to such solar neutrino scattering events and can probe large parts of the region of parameter space relevant for explaining (g − 2) µ . In particular, both experiments can probe the region where the preferred bands of (g − 2) µ and H 0 overlap. These are especially promising prospects in the case of LZ as this is a G2 experiment which will happen in the near future. In this context, we want to note that this interesting region can also be probed entirely in K → µ3ν at NA62 and in a future COHERENT run, while neutrino trident production at DUNE still has some sensitivity.
In the lower panel of Fig. 6 we show the corresponding projections for the direct detection experiments derived with the neutrino fluxes from a low metallicity Sun. Overall the sensitivities are slightly worse due the generally lower flux of 7 Be and 8 B neutrinos. Nevertheless, the general picture does not differ substantially from the high metallicity case. The most noticeable difference is that due to the weaker Borexino limits, a large part of the parameter space relevant for a simultaneous explanation of (g − 2) µ and H 0 is still viable.
The results of Fig. 6 clearly show that the optimal search strategy to cover the (g −2) µ solution in the U (1) Lµ−Lτ model with direct detection experiments is to look for coherent neutrino-nucleus scattering and push for large exposures. For this range of mediator masses, electron recoils are a subdominant contribution, mostly due to the related background. Also, lowering the energy threshold is only effective for probing small mediator masses with nuclear recoils (see the rightmost panel in Fig. 4). For this reason, xenon-based experiments as LZ outperform low-threshold detectors such as SuperCDMS.

V. CONCLUSIONS
Models of gauged U (1) Lµ−Lτ symmetry have received a lot of attention in recent years since their special flavour coupling structure can help to explain tensions in the observed value of the anomalous magnetic moment of the muon (g − 2) µ and of the Hubble parameter H 0 . Motivated by these intriguing prospects, in this article we have studied the potential sensitivity of solar neutrino scattering to probe these conjectured explanations within U (1) Lµ−Lτ . We have summarised our main results in Figs. 5 and 6.
• We have reevaluated previous results from solar neutrino-electron scattering at Borexino. We find that these have been previously underestimated (see Appendix C for a similar result for U (1) B−L ) and that they yield the most constraining current limit in the region where (g − 2) µ can be explained. In the scenario of a high metallicity Sun, which is currently favoured by experimental data, the updated limits exclude a substantial part of previously untested parameter space relevant for a simultaneous explanation of (g − 2) µ and H 0 . In the case of a low metallicity Sun, Borexino cannot exclude large parts of the simultaneous explanation of (g − 2) µ and H 0 , but still yields constraints which are competitive with previous leading bounds.
• We have analysed recent data of the liquid argon run of the COHERENT experiment and derived the corresponding limits on the parameter space of a U (1) Lµ−Lτ boson. While these limits are similar in reach to those derived from a previous caesium-iodide run they are slightly less constraining and cannot exclude any untested parameter space.
• In view of the recently reported excess of electron recoil events in XENON1T, we have shown that a U (1) Lµ−Lτ boson cannot account for the observed excess. Instead, we have used the corresponding data to derive a constraint on the parameter space of U (1) Lµ−Lτ . We have found that this does not yield a competitive bound compared with already existing limits.
• Furthermore, we have shown that upcoming and future dark matter direct detection experiments are a very promising means of probing the phenomenologically interesting region of parameter space where a simultaneous explanation of (g − 2) µ and H 0 is possible. The main search strategy to cover this region is to look for coherent neutrino nucleus scattering with increasing exposures, while electron recoils and lower thresholds are not as optimal. Thus, while SuperCDMS and DarkSide-20k will not be able to set competitive limits without major improvements of the detectors, LZ and DARWIN will be able to test a significant part of the parameter space relevant for (g − 2) µ . Considering that LZ is a second generation DM experiment and will happen in the near future these are very exciting prospects.
We eagerly await the upcoming results from the measurement of the muon anomalous magnetic moment at the Fermilab E989 experiment. Should there be a confirmation (5 σ discovery) of the excess over the SM value, our work suggests that direct detection experiments would play a leading role in identifying the new physics model by probing the U (1) Lµ−Lτ parameter space in complementarity with COHERENT, Dune and NA62.
Note: Upon completion of this work we became aware of a very recent article [133] complementary to our analysis that explored modifications of the CEνNS neutrino floor in a gauged U (1) Lµ−Lτ .
In order to maintain full generality, we consider a set of (chiral) fermions ψ i carrying charges Q i a and Q i b under two Abelian gauge groups U (1) a × U (1) b according to with D µ = ∂ µ + i g a Q i a X aµ + i g b Q i b X bµ . At the loop level the fermions ψ i will induce a vacuum polarisation term for the two bosons X a and X b via the diagram which is analogous to the vacuum polarisation of the photon in QED. Summing over all the fermion fields ψ i running in the loop, the amplitude of the vacuum polarisation diagram is given by In order to evaluate the integral over the loop momentum we follow standard textbook methods 16 and introduce Feynman parameters to bring the integral into the form Completing the square and shifting the loop momentum to = k + xq we will finally end up with the expression for the vacuum polarisation as (A5) The momentum integral is clearly divergent (∼ d 4 (1/ 2 )) and we have to regularise it. We will use dimensional regularisation in d = 4 − ε dimensions to compute the divergent integrals. The result of this calculation is given by where µ denotes an arbitrary mass scale and ∆ i = m 2 i − x(1 − x)q 2 . In an MS renormalisation scheme we will split off the divergent 2/ε piece and absorb it into a counterterm [139]. Extracting the Lorentz tensor structure the finite renormalised part of the vacuum polarisation is then obtained to be The effect of this vacuum polarisation term will manifest itself through an effective mixing of the two gauge bosons X a and X b via In order to see how this comes about we have to match the diagram of the vacuum polarisation with the one originating from the mixing, This matching allow us to identify the effective mixing parameter as

Appendix B: Matter-induced Solar Neutrino Oscillations
For a model with non-flavour-universal couplings to neutrinos, it is important to know the population of each neutrino eigenstate incident on a potential experiment. To derive constraints from the scattering of solar neutrinos, we must therefore understand the matter-induced neutrino oscillations which take place inside the Sun. Although solar neutrinos are produced through multiple processes, they are all generated in the electron neutrino eigenstate. This means that the fraction of the neutrino flux incident on our detector in the α flavour eigenstate will be given by the transition probability P (ν e → ν α ) [140].
So far, the only direct observations of solar neutrinos have been through interactions with electrons [97][98][99]. In the standard model, at the typical energies of solar neutrinos, it is valid to assume that the cross section for a muon-and tau-neutrino scattering with an electron is identical, as both proceed through the exchange of a Z-boson (the "neutral current" channel), whereas the scattering of electronneutrinos has the additional "charged current" contribution from the exchange of a W -boson. Under this assumption, the only transition probability that must be known is the electron-neutrino survival probability, P (ν e → ν e ). This quantity can be calculated with high precision in the two-neutrino approximation [140], and this has been the subject of many previous works [141][142][143][144]. However, in our model each flavour eigenstate has a different coupling, so it is important to know the relative abundance of each eigenstate.
As discussed in Section II.2, previous works have shown that any non-flavour diagonal NSIs can affect the matter potential and can therefore change the neutrino osciallation probabilities [91]. However, as per Eq. (30), assuming that solar matter is electrically neutral means that the NSIs couplings generated in our model with protons and electrons exactly cancel and no net effect is expected on solar neutrino oscillations.
Oscillation probabilities will also differ between the various solar neutrino fluxes, since they are produced in different regions of the Sun [144]. Fig. 7 shows the averaged three-flavour oscillation probabilities for all the relevant solar neutrino fluxes. These values were calculated using the nuSQuIDS software [145]. boson, the differential neutrino-electron scattering cross section including interference is found to be, where the couplings g α 1 and g α 2 are given by Eq. (27). While in the derivation of [46] these interference effects have been taken into account, important theoretical uncertainties on the SM and SSM prediction (including the uncertainty stemming from the solar metallicity problem) have been omitted. Combined with the consistent incorporation of the interference term in Eq. (C2) the inclusion of these uncertainties warrants for an updated analysis of the corresponding Borexino constraint on U (1) B−L as outlined in Section III.1.
We show the results of our improved analysis in Fig. 8. Existing constraints on the parameter space of a U (1) B−L boson are shown as the grey areas 17 . In the top panel we show the reevaluated Borexino  − 2)µ. The yellow dotted line is the upper limit derived form the recent XENON1T result [52]. Top: The red area corresponds to the Borexino exclusion limit in the scenario of a high metallicity (HZ) Sun. Bottom: As above, but with the Borexino limit calculated in the scenario of a low metallicity (LZ) Sun. See text in Appendix C for detailed explanations. limit from the 2011 data set [109] as the red dashed line for the case of a high metallicity Sun. The constraint derived from the more recent 2017 data set [108] is shown as the red solid line. Interpreting these two lines as an envelope of this true constraint, which would require the two results to be combined in a single chi-squared analysis with appropriate handling of any shared systematic uncertainties, we see that Borexino and NA64 [146] are setting the most stringent bound for moderate couplings in the range M A ≈ 20 − 80 MeV. For comparison we show the previous bound derived in [46] as the red dotted line.
We also show the results of our simplified analysis of the recent XENON1T result, discussed in detail in Section III.3.1. Although we have not performed the most complete analysis, we expect it to give an optimistic estimate of the limits, and since they do not cover any new regions of the parameter space we do not expect this result from XENON1T to supply competitive constraints on a U (1) B−L gauge boson within this mass region.
In the bottom panel of Fig. 8 we have evaluated the limits from Borexino and XENON1T in the case of a low metallicity Sun. In this case, the lower flux of 7 Be neutrinos leads to a weaker constraint from Borexino than we obtain in the HZ case, and we do not exclude any more of the parameter space beyond the limits already set by NA64 [146].
Finally, we expect qualitatively quite similar results of an updated analysis of the Borexino data for other anomaly-free groups gauging a combination of B and L e , like U (1) B−3Le , U (1) B−Le−2Lα or U (1) B−Lα−2Le with α = µ, τ .