On $(g-2)_\mu$ From Gauged $\mathrm{U}(1)_X$

We investigate an economical explanation for the $(g-2)_\mu$ anomaly with a neutral vector boson from a spontaneously broken $\mathrm{U}(1)_X$ gauge symmetry. The Standard Model fermion content is minimally extended by 3 right-handed neutrinos. Using a battery of complementary constraints, we perform a thorough investigation of the renormalizable, quark flavor-universal, vector-like $\mathrm{U}(1)_X$ models, allowing for arbitrary kinetic mixing. Out of 419 models with integer charges not greater than ten, only 7 models are viable solutions, describing a narrow region in model space. These are either $L_\mu-L_\tau$ or models with a ratio of electron to baryon number close to $-2$. The key complementary constraints are from the searches for nonstandard neutrino interactions. Furthermore, we comment on the severe challenges to chiral $\mathrm{U}(1)_X$ solutions and show the severe constraints on a particularly promising such candidate.

There are two immediate questions pertaining to the anomaly-free U(1) X gauge extensions of the SM: i ) Is U(1) Lµ−Lτ the only phenomenologically viable possibility that can explain the (g − 2) µ anomaly?
ii ) If there are alternative models, can these be experimentally disentangled from U(1) Lµ−Lτ ?
In this paper, we systematically explore the two questions, building on our previous work, Ref. [62]. The main result is that the only significant deviation of L µ − L τ allowed by data is gauge groups where L e −2B and the kinetic mixing with the photon approximately cancels the electron charge. This conclusion is rather nontrivial since there exists an extensive list of precise experimental probes constraining complementary combinations of 1 The issue of the SM prediction may not be completely settled. The lattice QCD determination of the hadronic vacuum polarization contribution to (g − 2)µ by the BMW group reduces the tension between the experimental world average and the SM (g − 2)µ prediction to less than two standard deviations [25]. It is expected that other lattice groups will be able to weigh in on this issue in the near-to medium-term future.
2 Both transitions are due to dipole moment operators, L eff ⊃ −e v¯ i L σ µν j R Fµν /(4πΛij) 2 + H.c., where v = 246 GeV is the electroweak vacuum expectation value. The NP contribution to (g − 2)µ is due to a flavor-conserving dipole, i = j = 2, whereas µ → eγ decay is from a flavor-changing dipole, i = 1, j = 2 or i = 2, j = 1. U(1) X gauge boson couplings, resulting in a limited number of phenomenologically viable possibilities.
The main working assumptions for our analysis are i) the SM is minimally extended by three right-handed neutrinos and a flavor-non-universal U(1) X gauge symmetry and ii) the theory is anomaly-free. The NP models also contain an SM singlet scalar field φ charged under U(1) X , whose vacuum expectation value (VEV) v φ breaks the U(1) X . The solutions to the (g − 2) µ anomaly typically require v φ around the electroweak scale if muons and φ carry O(1) U(1) X gauge charges. The details about the source of U(1) X symmetry breaking are not important for most of the phenomenology discussed in this paper, so our results do not change if a different SM-singlet condensate breaks U(1) X (or several condensates). Actually, viable neutrino masses and mixings typically require more than one scalar field.
Requiring the ratios between (non-zero) U(1) X charges for the chiral fermions to be at most ten gives roughly 21.5 million inequivalent integer charge assignments, up to flavor permutation [63]. We restrict our analysis to a subset of these, the 255 quark flavor-universal vector-like U(1) X charge assignments [62]. Taking into account the flavor permutations among charged SM fermions this gives 419 phenomenologically distinct U(1) X models. Here, we assumed that the sterile neutrinos are heavy enough so as not to affect the low-energy phenomenology, and, thus, different charge assignments for sterile neutrinos lead to the phenomenologically equivalent U(1) X model in this counting. In Section 5, we furthermore comment on the 21 chiral U(1) X charge assignments. The choice of quark flavor universality of U(1) X charges is phenomenologically well motivated, since flavor-non-universal U(1) X charges of quarks lead to dangerous flavor-changing neutral currents (FCNCs). As illustrated in Ref. [62] by a benchmark example, even for the second-safest option, a third-quark-family model with down-alignment, the CKM rotation induced up-sector FCNCs effectively rule out the parameter space preferred by (g − 2) µ . 3 While quark flavor-universal models avoid FCNCs, they are still severely constrained through a combination of other measurements: i ) neutrino trident constraints [51,52,65,66]; ii ) electroweak precision observables [67][68][69]; iii ) neutrino oscillation constraints on nonstandard neutrino interactions (NSI) [70,71]; iv ) measurements of coherent neutrino scattering on nuclei [72][73][74]; v ) the Borexino measurement of the cross section for the elastic scattering of 7 Be solar neutrinos on electrons [75,76] and other elastic neutrino-electron scattering experiments [77][78][79][80]; vi ) searches for new light resonances [59,81].
In Ref. [62], we studied the implications of these measurements for several selected benchmarks. In this manuscript, we go well beyond the initial analysis of Ref. [62] and assess the constraints for the complete set of 419 distinct vector-like U(1) X models, focusing on the parameter space region relevant for explaining (g − 2) µ . In several instances our phenomenological analysis applies to the complete class of quark flavor-universal vector-like charge assignment, even allowing for arbitrarily large charge assignments. We also comment on the phenomenology of chiral charge assignments, working out the details of theL µ−τ model in which the couplings of X µ to electrons are purely axial, thus, eliminating the very strict constraints on NSI from neutrino oscillations.
In the bulk of the paper, the analysis is kept as general as possible. In particular, we do not impose any requirements regarding possible connections with the b → s anomalies. This does not mean that such connections do not exist. On the contrary, lepton flavornon-universal U(1) X gauge symmetries can further support solutions of the B-physics anomalies that rely on tree-level exchanges of leptoquarks (LQ). In general, TeV-scale LQs tend to excessively violate the accidental symmetries of the SM-baryon and individual lepton number symmetries-all of which are exquisitely tested experimentally. Charging the LQ under a flavor-non-universal U(1) X gauge symmetry can reinstate the accidental symmetries while keeping the contributions to B-anomalies intact. Prominent examples of such mediators are the muoquarks [82][83][84][85][86], LQs charged under a U(1) X such that they carry global muon and baryon numbers. Both of these remain accidentally conserved at dimension 4, as they are in the SM. Out of 255+21 quark flavor universal charge assignments, 252+21 satisfy the muoquark criteria for the scalar weak triplet mediator [62].
The paper is organized as follows. In Section 2, we introduce the anomaly-free U(1) X models and the parameter space relevant for (g − 2) µ bounded by cosmology (m X 10 MeV) and perturbative unitarity (m X 1 TeV). Section 3 contains the discussion of the experimental constraints relevant for vector-like quark universal U(1) X models. The phenomenological implications of these constraints are presented in Section 4: In Section 4.1, we use the neutrino trident production, combined with the Z-pole constraints, to set a robust upper limit on the X boson mass (m X 4 GeV) for all renormalizable models introduced in Section 2. In Section 4.2, we perform a global analysis of experimental constraints and identify seven viable vector-like models that can explain the (g − 2) µ anomaly. Section 5 contains a brief discussion of chiral models and a detailed phenomenological analysis of the most promising example, theL µ−τ model. Possible connections with the B anomalies are discussed in Section 6, while Section 7 contains conclusions. The details on the calculation of neutrino oscillation bounds are given in Appendix A, while Appendix B contains details on the construction of the global χ 2 function used in the analysis.

Model framework
We start by reviewing the salient features of the gauged U(1) X models and how these can address the (g − 2) µ .

Field content and U(1) X charges
The models we consider have the SM matter content extended by three sterile neutrinos N i while the SM gauge group is enlarged by an additional U(1) X factor. After electroweak symmetry breaking (EWSB), the part of the Lagrangian relevant to the phenomenology at low energies (m X m Z ) is given by 4 where F µν and X µν are the electromagnetic and U(1) X field strength tensors, ε is the kinetic mixing parameter, J µ EM the electromagnetic current, and J µ X = f x f f γ µ f the current associated with U(1) X , where x f are the U(1) X charges for the chiral field f . The kinetic mixing term ε 2 F µν X µν can be removed by performing a non-unitary transformation of the Abelian gauge fields [87], after which the U(1) X Lagrangian is given by setting ε → 0, g X → g X = g X / √ 1 − ε 2 in Eq. (2.1) and shifting the charges according to in the expression for J µ X (Q EM f is the electric charge of the SM fermion f ).

Charge assignments
The U(1) X charges of quarks are assumed to be universal, such that the quark Yukawa interactions are given by the usual dimension-4 operators. Without loss of generality, we can set the U(1) X charge of the SM Higgs to zero, x H = 0. This leaves 276 inequivalent charge assignment, not counting flavor permutations, with integer charges for the SM fields in the range from -10 to 10, as listed in [82]. Shifting all the x f charges by a multiple of the hypercharge Y f gives physically equivalent models that would have x H = 0, see Appendix A.1 of Ref. [82]. The above quark flavor-universal U(1) X charge assignments fall into one of two categories. The 255 charge assignment in the vector category have vector-like U(1) X charges for both the quarks and the charged leptons, x L i = x E i , i = 1, 2, 3. The charged lepton masses are thus also generated via dimension-4 SM Yukawa interactions. Viable neutrino masses and mixings usually require additional U(1) X -breaking scalars (SM singlets), which lead to Majorana masses through N i N j φ ij interactions. We assume that the sterile neutrinos are heavy enough that they are not relevant for the low energy phenomenology. Taking into account the flavor permutations of the vector-like U(1) X charge assignments this leaves us with 419 phenomenologically distinct vector-like U(1) X models. 5 The possible charge assignments for the models in the vector category are given by [82,88] where {B, L e , L µ , L τ } are the usual values of baryon and lepton numbers for the SM fermion f , while L N i are the right-handed neutrino numbers (N i are not charged under 3) need to satisfy the Diophantine equations [63,89] (c e + c µ + c τ = i c N i and c 3 e + c 3 µ + c 3 τ = i c 3 N i ) giving a 4 parameter family of models valid beyond the restriction to charges less than 10. In particular, there are valid solutions for any set of c e , c µ , c τ : c N i = (−c e , −c µ , −c τ ). The parameterization in Eq. (2.3) will be exploited in Section 3.3 to asses NSI constraints on a large set of models. The charge assignments c N 1 = −c N 2 , c N 3 = c e,µ,τ = 0 (and permutations thereof) correspond to the dark photon type solutions to (g − 2) µ [90]. In this case, the X couplings to the SM fermions are exclusively due to kinetic mixing with the photon-typically radiatively induced from X interacting with some hidden sector particles also charged under the SM gauge group. This scenario has been ruled out as the solution of the (g − 2) µ anomaly, both when X decays visibly [91] or invisibly [92,93], with the possible exception of a combined decay [94].
There are additional 21 charge assignments for which there is no permutation such that x L i = x E i for every i = 1, 2, 3. These models constitute the chiral category and are listed in Section 2.2.2 of [82]. We examine their phenomenology and possible relevance for (g − 2) µ separately, in Section 5. In chiral models at least some of the charged lepton Yukawa couplings to the Higgs are forbidden at dimension four, leading to Yukawa matrices of rank less than three. Hence, some of the charged lepton masses are generated through higherdimension operators. We assume that these operators are either generated by integrating out vector-like fermions (which do not change the conditions for anomaly cancellation), or by integrating out heavy scalars. We discuss this in more detail for theL µ−τ model, including the phenomenological constraints, in Section 5.

Mater unification
Interestingly, some of the above U(1) X extensions of the SM can be unified into a larger semi-simple gauge group at higher energies. As an example consider the vector category charge assignments with x i = x N i . Ref. [85] showed that in that case the U(1) Y × U(1) X can be unified into SU(2) R × U(1) B−L × SU(3) lep . Starting from this semi-simple group the unification can proceed further (see Figure 1 of Ref. [85]).
2.4 Explanation of (g − 2) µ A massive vector X µ coupling to muons with interaction gives rise to a one-loop radiative correction to the anomalous magnetic moment of the muon: where r µ = m µ /m X . For the full 1-loop expressions see, e.g., Refs. [62,95,96]. Correction (2.5) is of the right sign to explain the difference between the measured and predicted anomalous magnetic moment of the muon, ∆a µ = a avg µ − a SM µ = (251 ± 59) × 10 −11 [2], if X µ couples mainly vectorially to muons.
Requiring that ∆a µ is explained by a massive vector X µ at one-loop order translates to an upper bound on the U(1) X breaking VEV , which is saturated for q A = 0. Perturbative unitarity, g X q V ≤ √ 4π, then implies an upper bound on the X µ mass, m X 1.0 TeV (see Ref. [97] for a more detailed discussion). This is a rather restrictive bound, the origin of which can be traced back to the fact that in the one-loop X µ contribution to (g − 2) µ , the required chirality flip necessarily occurs on the muon leg and, thus, is suppressed by the small muon mass.

Cosmology
A robust lower limit on the gauge boson mass m X follows from the agreement between observations of primordial light element abundances and the predictions within standard cosmology Big Bang Nucleosynthesis (BBN). The gauge coupling g X required to explain the (g − 2) µ anomaly is large enough that the gauge boson X µ efficiently thermalises with the SM plasma in the early universe, i.e., for temperatures T > m X . A light X µ contributes with additional relativistic degrees of freedom during BBN, changing the predictions for light element abundances. This is avoided for m X 10 MeV, in which case X µ decays quickly, before the onset of BBN. The precise bound on m X was derived for the U(1) Lµ−Lτ model in Ref. [98] (see also [60]), and holds approximately also for all the other U(1) X models considered here.

Experimental constraints
As a next step in the analysis, we discuss all the various constraints on the U(1) X models. They include EW precision tests, a variety of neutrino interactions, and finally resonance searches.

Electroweak precision tests
After EWSB the U(1) X gauge boson X µ mixes both with the photon, A µ , and the Z boson. The mixing with photon is important for low energy constraints, while the mixing with the Z is severely constrained by the electroweak precision tests.

Electroweak symmetry breaking
In the models we consider, the source of the EWSB is the same as in the SM-the VEV of the SM Higgs. Above the electroweak scale the kinetic mixing Lagrangian between X µ and the photon, Eq. (2.1), is replaced by the kinetic kinetic mixing between X µ and B µ , parametrized by the parameter ε Y : where D µ H = (∂ µ + ig X x H X 0,µ )H. In writing the above Lagrangian we remain agnostic about the origin of the X boson mass. The mass term m X 0 could be due to a VEV of a single SM singlet scalar or due to a more complicated condensate in the SM singlet sector. The U(1) X charge of the SM Higgs, x H , can be absorbed in the other parameters by performing the shift where ξ = 2x H g X /g 1 . This shift changes the charges x f of the matter fields, f , by with Y f the hypercharge of fermion f . Without loss of generality, we can, therefore, take x H = 0, which is what we will assume from now on. Expressing B µ and W 3 µ in Eq. (3.1) in terms of the photon, A µ , and the Z boson field, Z µ , gives the kinetic mixing Lagrangian after EWSB: where m Z 0 = g Z H / √ 2 is the SM Z mass. In Eq. (3.4) we do not write down the terms involving dynamical Higgs field. We also used the shorthand notation s w = sin θ w , c w = cos θ w , and t w = tan θ w , where θ w is the weak mixing angle.
A simultaneous diagonalization of the kinetic and mass terms is achieved by a combination of a non-unitary field redefinition and a rotation, . The mixing angle θ between Z and X is given by where with m X the physical X boson mass. The diagonalization of the gauge fields changes their effective currents. From Eq. (3.5), we find that where g X = c θ ∆ ε g X is the physical X gauge coupling.

Z mass constraint
The kinetic mixing between X µ and B µ introduces a mass mixing between Z and X gauge bosons, resulting in a shift of the physical Z mass m Z and corrections to Z couplings to SM fermions. While a global fit to the electroweak precision data would be required to capture all the resulting experimental constraints on X µ -B µ mixing, it suffices for our purposes to consider just the T parameter (marginalized over other electroweak observables), mainly because the measurements of SM fermion couplings to Z have larger relative errors, see also discussion in Ref. [99], as well as Refs. [100][101][102][103], and Section 3.1 below. We find that where m Z,SM is the SM prediction for the Z mass (from the W mass including radiative corrections), whereas m Z,obs is the measured Z mass so that in the SM ρ 0 = 1. Since we are interested in the NP constraints, we can use the tree-level relations, giving the second (approximate) equality. ρ 0 is related to the oblique parameter T through ρ where experimentally from the electroweak fit T = 0.03 ± 0.12, when the S, U parameters are allowed to float freely [67] (see also [104]). The Z mass shift in our model results in (3.12) In the limit of small X mass and small kinetic mixing parameter, we have This provides a relatively weak bound for the region of light X masses (m X < 2m µ ) well above the typical expectation ε Y ∼ O(eg X /16π 2 ) for radiatively induced kinetic mixing. Nevertheless, the T parameter bound in Eq. (3.14) is phenomenologically quite important, since it does not allow for arbitrarily large kinetic mixings. Combined with constraints from neutrino trident production, it translates to a model-independent requirement that X needs to be lighter than a few GeV (see Section 3.2 for details). For X masses comparable to m Z , the bound in Eq. (3.14) becomes more stringent and then progressively weaker for m X m Z . Throughout this region, the mixing angle θ between X and Z is given by the upper expression in Eq. (3.6), except for the very small region where X and Z are almost mass degenerate, such that |s wεY | > |1 − r 2 X |. In this mass degenerate regime, the constraint on the T parameter leads to As anticipated, in this region, the constraint on |ε Y | is significantly more stringent than it is for the light X mass limit, Eq. (3.14). From the T bound, it follows that this case is only relevant for |1 − r 2 X | < 2.2 × 10 −3 , i.e., when the X and Z are degenerate to within 0.2%.

Lepton universality in Z decays
The mixing between X and Z also results in non-universal Z boson couplings to the SM leptons, Eq. (3.9), which were measured to high precision at LEP [69]. We expect the strongest constraint to come from the lepton flavor universality ratio for first two generations of leptons [69], Using Eq. (3.9), we obtain for the leading new physics contributions to flavor universality ratio. This bound is more model-dependent than the T bound, since it depends on the U(1) X charges of the leptons that change from one U(1) X model to another. To be completely consistent, one would in principle have to perform a global electroweak fit for every U(1) X model. However, we estimate the typical non-universal effects to be sub-leading and, thus, their inclusion would only marginally improve on the constraint obtained from the T parameter in Section 3.1.
Effective X current After EWSB the effective couplings of X boson (the mass eigenstate mostly composed of X 0 ) are encoded in the effective current J X,eff , Eq. (3.10). In the limit |s wεY | < |1 − r 2 X | the current takes a simpler form: For small masses, m X m Z , the last term is power suppressed and we find agreement with the low-energy description of X mixing with the photon, given in Section 2.

Neutrino trident production
The ∆a µ anomaly points to a vector boson in the mass range 10 MeV m X 1 TeV. The viable mass window is set by cosmology (lower limit, see Section 2.5) and perturbative unitarity (upper limit, see Section 2.4). An efficient complementary constraint, which cuts significantly into this parameter range, is due to limits on nonstandard neutrino trident production, i.e., the scattering of a muon neutrino on a nucleus, producing a pair of charged muons, ν µ N → ν µ N µ + µ − . In combination with the constraints from the electroweak T parameter, Section 3.1, it limits the X mass to m X 4 GeV, as we show below. Neutrino induced production of a µ + µ − pair in the Coulomb field of a heavy nucleus is a rare electroweak process in the SM. In U(1) X models, there is an additional tree-level contribution from the diagram with an X µ gauge boson exchanged between ν µ and µ legs. The strongest bound on this contribution is due to the CCFR experiment [65], which reported the measurement for a trident cross section normalized to the SM prediction. 6 This imposes constraints on the vector and axial vector couplings of X µ to muons. We derive the resulting bounds on the gauge coupling g X as a function of m X for the U(1) X models using the public code of Ref. [52] (further details can be found in Ref. [62]). In the EFT region, m X 1 GeV, this bound is approximated by [52] σ CCFR σ SM are the normalized NP coefficients of the effective 4-fermion interactions between muon neutrinos and vector and axial muons, respectively. G F is the Fermi constant, controlling the overall strength of the SM contribution, while q V,A and q ν are the effective U(1) X vector and axial charges of muons and the charge of muon neutrinos, respectively (cf. Eq. (2.4)):

Neutrino oscillations and coherent elastic neutrino-nucleus scattering
Beyond the trident production, discussed in Section 3.2, the U(1) X solutions to the (g − 2) µ anomaly also lead to two other types of nonstandard neutrino interactions (NSI): the modified matter effects in neutrino oscillations and the additional contributions to coherent elastic neutrino-nucleus scattering. The matter effects in neutrino oscillations are given by the forward scattering amplitude, i.e., at zero momentum transfer. Accordingly, the NSI contributions are well described by the EFT Lagrangian (f = {e, p, n}) [71,[106][107][108], even for m X well below the mass window of interest, m X 10 MeV. In our setup, the NSI are generated at tree level by integrating out the gauge field X µ . We use the results of a global EFT fit to the neutrino oscillations data [70,71,109] to set constraints on X µ couplings to u and d quarks and to electrons (cf. App. A). The constraints are numerically important for the (g − 2) µ compatible parameter space. For instance, the upper limit on g X /m X from neutrino oscillations rules out U(1) B−3Lµ as a possible solution to the (g − 2) µ anomaly, despite relatively small couplings to quarks, x q /x µ 1 (see Fig. 3 in Ref. [62]). The bounds on NSI from neutrino oscillations, Eq. (3.24), are sensitive to the average U(1) X charge of the material the neutrinos propagate through, either the Earth's mantel and core or the sun. Among other implications, this means that the value of the kinetic mixing ε has no effect on the neutrino oscillation bounds, since the matter is electrically neutral. The bounds from neutrino oscillations are relaxed for a particular set of U(1) X models, such as B − 2L e − L τ and B − 2L e − L µ , for which the average U(1) X charge of normal matter is almost zero (nuclei have, on average, roughly as many protons as neutrons, i.e., B 2L e for normal matter).
A complementary set of constraints on NSI is due to coherent elastic neutrino-nucleus scattering [72,73], which was observed by the COHERENT experiment [74]. In this case, the EFT description in Eq. (3.24) is not valid for the full range of X masses of interest to our analysis, 10 MeV m X 4 GeV. Following Ref. [62], we instead keep X µ as a dynamical field when setting the bounds on g X , using the likelihood from Ref. [110].

Elastic neutrino-electron scattering
The Borexino experiment measured the cross section for elastic scattering of solar neutrinos on electrons [75,76], which constrains possible new X µ mediated interactions between neutrinos and electrons. The strongest bound on light vector boson interactions is obtained from 7 Be neutrinos, which have an energy of 862 keV. The direct X coupling to electrons can significantly change the scattering cross section of solar neutrinos, especially due to many U(1) X models having large couplings to muon and tau neutrinos. The ν µ and ν τ neutrino flavors are present in the solar flux on Earth, since the initial ν e neutrinos oscillate during propagation from the Sun. We follow the analysis of Ref. [52] to place bounds on g X while requiring that the scattering cross section remains within 2σ of the measurement.
The reactor experiments TEXONO [77] and GEMMA [78] measured the related cross section for elastic scattering of electron anti-neutrinos on electrons, while the high-energy beam experiment CHARM-II at CERN measured the cross sections for elastic ν µ e − and ν µ e − scattering [79,80]. These bounds can be as relevant as Borexino and will be discussed in the context of the chiralL µ−τ model in Section 5, where the electron coupling is purely axial avoiding bounds from NSI oscillations.

Resonance searches
Several intensity-frontier collider experiments [59,81] have directly searched for a production of a vector resonance X in the mass range of interest, 10 MeV m X 4 GeV.
In a fixed target experiment such as NA64 [111], the vector boson is produced via bremsstrahlung process eN → eN X, where N is a nucleus. The most relevant decay mode is X → invisible, which is typically dominant below the dimuon threshold. The events are reconstructed from the missing energy measurements. The NA64 search [111] constrains the X couplings to electrons, which enter the prediction for the X production rates. Future NA64µ [112] and M 3 [113] experiments will feature muons in the incoming beam and will have the potential to cover the entire parameter space of the L µ − L τ model, see Fig. 2 in Ref. [62]. The fixed target experiments are effective only for m X 1 GeV.
Another important set of direct searches was performed at B-factories and probed X masses up to 10 GeV. The BaBar search for e + e − → µ + µ − X in the 4µ final state [114] (labeled in figures as BaBar 2016) sets fairly stringent constraints above the dimuon threshold. BaBar also searched for a radiative return process e + e − → γX with X decaying to e + e − or µ + µ − [115] (BaBar 2014 and LHCb) or to invisible [93] (BaBar 2017). The LHC searches extend the exclusion to even larger X masses, such as the LHCb search for X → µ + µ − [116,117] and the CMS search for Z → µ + µ − X → 4µ [118].
In all the aforementioned searches, the X decays are prompt in the parameter space relevant for (g − 2) µ . We use the DarkCast code [81], which comes with the compilation of relevant bounds, to set limits on the gauge coupling g X as a function of the mass m X . Crucially, the above bounds are model dependent; for instance, the constraints from dimuon resonance searches could be removed by introducing additional invisible X decays to a light dark sector.

Phenomenological implications
We explore next the phenomenological implications of experimental constraints on minimal anomaly-free U(1) X models as candidates for explaining the (g −2) µ anomaly. In Section 4.1, we show that a combination of trident and electroweak precision tests, assuming the SM is the only source of EWSB, i.e., that U(1) X is broken by SM singlet scalar(s), leads to the upper bound m X 4 GeV. In Section 4.2, we include other experimental constraints that are particularly relevant for this low X mass region and perform a global analysis of the complete set of 419 phenomenologically distinct vector-like U(1) X models with charges up to 10. The discussion of chiral models is relegated to Section 5.

Generic upper bound on m X
Ref. [56] demonstrated that the trident production sets stringent constraints on the viable parameter space of the (g − 2) µ solution in the U(1) Lµ−Lτ model, limiting The above constraint was shown in Ref. [62] to hold for all the U(1) X models in the limit of vanishingly small kinetic mixing (ε = 0), where the muon couplings to X µ are due to muon U(1) X charges. This was achieved by marginalizing the trident constraint over all values of x L 2 and x E 2 -the U(1) X charges of the second generation left-handed lepton doublet and right-handed lepton singlet, respectively-with the results of the analysis shown in Fig. 1 of Ref. [62]. In the ε = 0 limit, the vector and axial couplings of X µ to muons are directly correlated to the U(1) X charges of left-and right-handed muons, with the vector charge given by (x L 2 + x E 2 )/2 and the axial charge by (x L 2 − x E 2 )/2. This is the reason why the trident bound so efficiently constrains the possible U(1) X solutions of (g −2) µ . Since (g −2) µ requires near-vectorial couplings, that is, comparable x L 2 and x E 2 , a considerable x L 2 is required. 7 This, in turn, implies similar couplings of X µ to muons and muon neutrinos, the upper components of the L 2 doublet, and, thus, a sizable neutrino trident production ν µ N → ν µ N µ + µ − . When ε = 0, the direct correlation between the NP shift in (g − 2) µ and the trident production no longer applies. The kinetic mixing between X µ and the photon in the lowenergy Lagrangian (2.1) shifts the coupling of upper and lower components of the isospin doublets independently and decorrelates muon and neutrino couplings to X µ . An instructive example is the B − 2L e − L τ model solution to (g − 2) µ , in which the muon couplings to X µ are generated entirely through its mixing with the photon, while there are no couplings to muon neutrinos and, thus, essentially no trident bounds on X µ . 8 Alternatively, for m X comparable to m Z , one might imagine a scenario where the J X and J Z contributions 7 In the renormalizable, anomaly-free U(1)X vector models, the Xµ couplings to charged leptons are diagonal in the mass basis and do not lead to cLFV. The new physics contribution to (g − 2)µ is due to a muon and an Xµ running in the loop. The bottom-up phenomenological vector model of Ref. [119], in which the new physics contribution to (g − 2)µ is due to a tau and an Xµ running in the loop, is difficult to UV complete in a gauged U(1)X . Such a solution requires, for instance, an extended scalar sector beyond our minimal setup (see, e.g., Ref. [120]). 8 The typical flavor composition of the initial neutrino flux in these experiments is dominantly νµ, whereas the next-largest, νe, constitutes less than a percent of the neutrino flux (see, e.g., Ref. [121]). Suppressing the X coupling to νµ is sufficient to make the CCFR bound irrelevant in the parameter region that leads to the solution of (g − 2)µ. Figure 1. Minimum χ 2 values in a model with a new X µ vector boson coupling to second generation leptons with charges x L2 and x E2 and a kinetic mixing with the hypercharge gauge boson ε Y from a fit to experimental data taking into account (g − 2) µ , CCFR, and the EW T parameter (see the text for details). The orange curve is the global best fit χ 2 , the blue is the best fit in the absence of muon charges (dark photon), the red is the best fit for vanishing kinetic mixing, and the green line is the SM χ 2 . In the shaded grey region, (g − 2) µ cannot be explained by any values of x L2 , x E2 , and ε Y while at the same time satisfying the constraints from CCFR and the electroweak T parameter.
to J X,eff (3.18) conspire to a similar effect. To account for these possibilities, we include EWPT as a complementary constraint to set a model-independent upper bound on m X for viable (g − 2) µ solutions.
Given m X , the X boson couplings to second generation leptons (3.18) are determined through a known combination of parameters g X x L 2 , g X x E 2 , and ε Y . These determine the U(1) X shifts in T (3.12), the neutrino trident cross section (3.21), and (g − 2) µ (2.5). We use these observables to construct a combined χ 2 . For each m X in 1 GeV − 300 GeV, we then find the χ 2 minimum by varying g X x L 2 , g X x E 2 , ε Y , shown as the orange line in Fig. 1, to be compared with the SM value in green. The case of no kinetic mixing, ε Y = 0, is shown with the red line, while the blue line shows the case where the couplings to X µ are entirely due to kinetic mixing.
We observe that at low m X masses, the best fit is similar to the solution where the X µ couplings are exclusively generated by the kinetic mixing. However, as the m X mass increases there is a growing tension between a large value of ε Y , which would minimize the size of the X µ couplings to neutrinos and the effect of the trident constraints, and a small value as preferred by EWPT constraints. The net effect is that a combination of trident production and EWPT constraints limit the X boson solutions to (g − 2) µ to be rather light: It is worth reiterating that the only assumptions entering this bound are that the U(1) X model is renormalizable. Higher-dimensional operators could be added to the Lagrangian (3.1) to modify the bound; however, this goes beyond the minimality assumed in this paper.

Global constraints on anomaly-free vector-like U(1) X models
Having set an overall upper bound on m X , we perform a scan over all vector-like models in order to determine those that can account for the (g − 2) µ discrepancy and at the same time satisfy all the constraints discussed in Section 3. We assume that the sterile neutrinos are heavy enough so as not to affect the low energy phenomenology. The U(1) X charges in the vector-like models then constitute a 3-parameter family of possible charge assignments, with the charges x f for each SM fermion given by with B, L e,µ,τ the baryon and individual lepton numbers and c e,µ,τ continuous parameters.
For the 419 vector-like models with |x f | ≤ 10, these take specific rational values that can be found in Ref. [62]. An overall factor can be absorbed into the gauge coupling g X , which reduces the number of independent physical parameters to two. Since it is well-known that the L µ − L τ model can explain (g − 2) µ while passing all constraints, we express the two-dimensional model space in Eq. (4.3) in a way that makes it apparent how closely a given model resembles the L µ − L τ charge assignments, The polar angle α determines the ratio of baryon and electron charges, while the radial parameter R gives the "closeness" to L µ − L τ , with R → ∞ corresponding to the exact L µ − L τ limit. Note that L e and B, with coefficients proportional to sin(α) and cos(α), respectively, enter the parameterization in a combination with L µ in such a way that models given by Eq. (4.4) are anomaly-free for arbitrary values of α and R. We use the parametrization (4.4) to systematically treat the experimental constraints on the vector-like models. Since only low X mass explanations of (g − 2) µ are still possible, cf. Eq. (4.2), we can focus on the low-energy Lagrangian (2.1). For each vector-like model, the physics is determined by the X mass m X , the gauge coupling g X , and the kinetic mixing parameter ε. To reduce the complexity of the analysis, we identify m X 200 MeV as the mass with the best possibility of explaining (g − 2) µ : for masses above the di-muon threshold the resonance searches at LHCb and Babar provide very strong constraints. All the remaining important bounds, in particular the bounds on NSI from neutrino oscillations, Borexino, NA64, and COHERENT, limit the X mass from below, and thus taking m X close to the di-muon threshold minimizes their importance.
To narrow down the possible candidates for explaining (g − 2) µ , we scan over α and R in Eq. (4.4). For each point in the (α, R) model space, we construct a ∆χ 2 function that takes into account the bounds from ∆a µ , Borexino, and NSI osc. and minimize it with respect to (g X , ε) at m X = 200 MeV. As an additional constraint, we require (g X , ε) to be in the range allowed by NA64 and BaBar2014 searches, as implemented in DarkCast. We deem a model to be excluded as a viable explanation of the (g − 2) µ anomaly, if the minimum ∆χ 2 value is above 4. This is to be compared with the SM value for ∆χ 2 , which is ∆χ 2 SM = 18.3), i.e., we discard models that are in tension with either (g − 2) µ or any individual constraint by more than 2σ. Further details on the construction of ∆χ 2 are given in Appendix B.
In the left panel of Fig. 2, we show the entire excluded region in the model space. For illustration, we mark all the models with integer charges |x f | ≤ 10 with small blue crosses and the three benchmark models B − 3L e , B − 2L e − L τ , and B − 2Le − L µ by orange, purple, and green tripods. As expected, for large values of the radial parameter, R 50, all the models are viable, irrespective of the value of α. All such models closely resembles the L µ − L τ model with only small deviations.
For smaller R values the viability of the models is strongly α dependent, i.e., dependent on the baryon-to-electron charge ratio. For small values of R, the viable models are restricted to a narrow region around tan(α) = −2/3 (depicted as the black dashed line in Fig. 2 left) so that the U(1) X electron charge is roughly −2 times the baryon charge. For such charge assignments, normal planetary matter, which has on average about as many protons as neutrons, is almost neutral under U(1) X . As discussed in Section 3.3, this results in the reduced relevance of neutrino oscillation constraints. By contrast, the U(1) X models in which normal matter carries a considerable U(1) X charge are ruled out as an explanation of (g − 2) µ precisely because of the excessive contributions to the NSI. There are only two ways for a model to both explain (g − 2) µ and predict the atoms of the most common elements to be essentially U(1) X -neutral. Either tan(α) ≈ −2/3 so the electron and nucleon charges approximately cancel, or R 1 so the electron and nucleon charges are small compared to the muon charge. The only model in which all atoms are exactly U(1) X -neutral is the L µ − L τ model, i.e., R → ∞. For models that differ considerably from L µ − L τ , i.e., for R 1, the cancellation of U(1) X charges in normal matter is the only way to avoid the stringent NSI bounds.
In the right panel of Fig. 2, we show the region around tan(α) ≈ −2/3 for R < 1. Again, we mark the models with integer charges |x f | ≤ 10 with small blue crosses and the three sample models B − 3L e , B − 2Le − L τ , and B − 2Le − L µ with orange, purple, and green tripods, respectively. In the lower half of the plot α is shifted by π compared to the upper half of the plot, and R is always positive. 9 We observe that for small values of R, the only models that can pass the neutrino oscillation bounds either lie in a very narrow region close to tan(α) = −2/3, where nucleon and electron U(1) X charges cancel inside atoms, or at around R ≈ 0.1 in the region −1 tan(α) −2/3.
The constraints from neutrino oscillations on NSI already exclude most of the parameter space for R 10. Another potentially relevant constraint in this region of model space is due to coherent elastic neutrino-nucleus scattering as measured by the COHERENT experiment. We impose this bound on all the models with integer charges |x f | ≤ 10 that are not already excluded by the neutrino oscillation data. (Since deriving the COHERENT bound is computationally expensive this constraints was not include in Fig. 2.) All the models that are not yet excluded before applying the COHERENT constraint appear in the right panel of Fig. 2 as small blue crosses on white background, apart from the L µ − L τ model, which is not shown, since it corresponds to R → ∞ limit. For all these models, we add the contribution from COHERENT to the previous ∆χ 2 function and minimize the newly obtained ∆χ 2 function at m X = 200 MeV with respect to (g X , ε). The charges and best-fit values of the models with ∆χ 2 < 4 are shown in Table 1. This table shows models belonging to three distinct classes, which we discuss separately next.
The L µ − L τ model The L µ − L τ model has the lowest ∆χ 2 value. It is able to explain (g − 2) µ comfortably without causing tensions with other experimental constraints. The best-fit point is found for vanishing kinetic mixing ε, however even for kinetic mixing values that one would typically obtain if this were generated at one loop, ε ∼ g X e/(16π 2 ), the L µ − L τ model can still explain (g − 2) µ without being excluded by other constraints. In Fig. 3 we show an example where the kinetic mixing is assumed to be induced by the muon and tau running in the loop.
The class 3B − 6L e − 3L τ + N(L µ − L τ ) with N ∈ {3, 4, 5, 6, 7} These models satisfy tan(α) = −2/3, i.e., they lie on the dashed vertical line in Fig. 2 and are all within the narrow region of model space where the U(1) X charges of electrons and nucleons in atoms approximately cancel. In Fig. 2, only the models with N < −2 were excluded. The addition of the COHERENT constraint now leads to the exclusion of all models with N < 3. As shown in Table 1, larger values of N generically correspond to lower ∆χ 2 values, that is, they can satisfy the constraints more easily. This is visualized in the left panel of Fig. 4, where we show the ∆χ 2 values for N ≥ 0 with the grey shaded region corresponding to ∆χ 2 > 4. The reason for the decrease in ∆χ 2 as N increases is that the muon charge increases as well and thus a smaller gauge coupling is requires to , 1, 2, 3, 4, 5, 6, 7}. The ∆χ 2 function takes into account the measurement of (g − 2) µ and the bounds from NSI osc., COHERENT, Borexino, NA64, and BaBar while the X mass was set to 200 MeV. The grey shaded region corresponds to ∆χ 2 > 4. The five models with ∆χ 2 < 4 are also listed in Table 1. explain (g − 2) µ . This can be seen in the sixth column of Table 1 and in the right panel of Fig. 4. Since an increase in N does not affect electron and baryon charges, the smaller values of the gauge coupling lead to weaker constraints from NSI. Ultimately, the models become more similar to the L µ − L τ model as N increases.
At the best-fit point, not only the ∆χ 2 value and the gauge coupling g X but also the kinetic mixing ε is correlated with N , which can be understood as follows: The COHERENT bound is due to a measurement of the coherent elastic neutrino-nucleus scattering in cesium iodide (CsI). Accordingly, the COHERENT bound becomes less constraining if the U(1) X charges of protons and neutrons in CsI approximately cancel. Cesium and iodine have neutron-to-proton ratios 78/55 and 74/53, respectively, such that the average neutron-toproton ratio in CsI is 1.4 to a very good approximation. In the 3B − 6L e − 3L τ + N (L µ − L τ ) class of models, the effective U(1) X charge of neutron and proton (taking into account the kinetic mixing ε) is 3 and 3 + ε e/g X , respectively. CsI therefore becomes approximately U(1) X -neutral for The Borexino bound, on the other hand, is due to a measurement of elastic neutrino-electron scattering and vanishes if the effective U(1) X charge of the electron is zero, i.e., when (4.6) Figure 5. Parameter space of viable vectorlike U(1) X models that can explain (g − 2) µ (the parameter space for L µ − L τ is shown in Fig. 3). The central value of ∆a µ is indicated by a dashed black line (the 2σ band by solid black lines). The exclusions by different measurements are color coded as given in the legend.
For small N , the relatively large gauge coupling needed to explain (g − 2) µ makes the Borexino bound very sensitive to deviations from Eq. (4.6). 10 As N increases and the gauge coupling decreases, larger deviations from Eq. (4.6) are allowed by Borexino and the best fit value of ε e/g X moves closer to Eq. (4.5) in order to better satisfy the COHERENT bound. This also means that the COHERENT bound is considerably more constraining at the best fit point for smaller N , which in turn leads to a larger ∆χ 2 value. The bounds on 3B − 6L e − 3L τ + N (L µ − L τ ) with N ∈ {3, 4, 5, 6, 7} are shown in Fig. 5 in panels a) through e), in the plane of the X mass m X and the gauge coupling g X , while for each plot the kinetic mixing ε is set to its best-fit value shown in Table 1. In these plots, we can observe the interplay of Borexino and COHERENT constraints described above. For small N , e.g., for N = 3 shown in panel a), g X has to be relatively large in order to explain (g − 2) µ . Consequently, the Borexino bound (shown in dark purple) requires ε to be very close to Eq. (4.6). This in turn leads to a relatively strong COHERENT bound (shown in red), which then excludes part of the parameter region preferred by (g − 2) µ in the m X mass range between 100 MeV and 200 MeV. Comparing this to the models with increasing N in panels b) through e), one observes that the size of g X preferred by (g − 2) µ decreases due to the increasing muon charge, the Borexino bound allows for smaller ε, and the COHERENT bound becomes less and less important until it becomes weaker than the bound from the neutrino oscillations for N = 7, cf. panel e).
The plots in Fig. 5 also demonstrate that for m X above the dimuon threshold, resonance searches, in particular from BaBar and LHCb, prevent an explanation of (g − 2) µ , whereas for lower values of m X , neutrino oscillation data provides the strongest bound. This justifies the choice of m X = 200 MeV used in Figs. 2 and 4 and Table 1.
It is interesting that the CCFR bound from neutrino trident production (shown in orange) becomes less constraining for smaller values of N and even completely disappears for N = 0. The reason is that the effective U(1) X couplings of the charged muon and muon neutrino are N − ε e/g X and N , respectively. For sizable ε and decreasing N , the region in parameters space where (g − 2) µ can be explained is, therefore, less and less excluded by the CCFR bound. In the large N limit, on the other hand, the model is closer and closer to the L µ − L τ model and the CCFR constraint becomes increasingly important.
The class 3B − 7L e − 2L τ + N(L µ − L τ ) with N = −1 These models satisfy tan(α) = −7/9 so that the cancellation of U(1) X charges inside atoms is less efficient than in the models with tan(α) = −2/3. Nevertheless, there is a window roughly around N = −1 where the bound from neutrino oscillations is relatively weak and an explanation of (g − 2) µ is possible (cf. right panel of Fig. 2). The virtue of this class of models is that both the Borexino and the COHERENT bounds can be comfortably satisfied. This is apparent from the constraints displayed in Fig. 5, panel f), which show that for m X below the dimuon threshold the neutrino oscillations provide the most relevant bound in the (g − 2) µ band, whereas COHERENT and Borexino are far less important. While the COHERENT bound vanishes for ε e/g X ≈ −7.2, as in Eq. (4.5), the effective U(1) X charge of the electron in the present class of models is −7 − ε e/g X such that the Borexino bound disappears for ε e/g X = −7. Consequently, values of ε e/g X ≈ −7 lead to very weak bounds for both Borexino and COHERENT. The best fit values for ε e/g X are indeed found to be very close to −7 by the global scan (see Table 1). The only model in the class 3B − 7L e − 2L τ + N (L µ − L τ ) that enters the list of viable models in Table 1 is the one with N = −1, but it actually reaches the second lowest ∆χ 2 of all the models in this list. This model is also notable for particularly weak bounds from neutrino trident production.

Chiral models
There are 21 chiral charge assignments in the charge range [−10, 10] (not counting permutations) in which lepton masses cannot be fully realised at the renormalizable level [82]. Upon careful inspection of the charge assignments in the entire class of the chiral U(1) X models, we focus on one of the best performing models, the modelL µ−τ . Quite generally, for these types of models the constraints in the parameter region relevant for (g − 2) µ will be minimized, if the U(1) X charges of the quarks vanish, and the electron coupling is purely axial, in order to avoid the NSI bounds, while the vector-like muon charge is as large as possible. TheL µ−τ model satisfies all of the above requirements. Even so, theL µ−τ model is excluded as the explanation of the (g −2) µ anomaly by a set of complementary constraints, as we show below, and we expect the same to be true for the other chiral models.
In theL µ−τ model, the U(1) X charges of the fermions are given by [62] ( where the first line is for the left-handed leptons, the second and the third lines are for the right-handed charged leptons and neutrinos, respectively, while the quarks do not carry a U(1) X charge. Due to the U(1) X charge assignment, the charged lepton Yukawas are forbidden at the renormalizable level and only arise once the U(1) X gauge symmetry is spontaneously broken. For concreteness let us consider the minimal possibility-that it is due to a SM-singlet scalar φ with U(1) X charge x φ = 1 that develops a VEV φ = 0. The diagonal entries of the charged lepton Yukawa matrix are populated by the dimension-5 operatorsL i Hφ † E i for muons and taus, i = 2, 3, and by the dimension-6 operatorL 1 Hφ 2 E 1 for the electrons. This is consistent with the smallness of the charged lepton masses and with the hierarchy between electron and muon and tau Yukawas. Such higher-dimension operators can be generated, for instance, by integrating out a set of heavy vector-like leptons at tree level with masses M φ well above the EW scale. The muon and tau Yukawas are then suppressed by φ /M and the electron Yukawa by ( φ /M ) 2 .
The off-diagonal terms in the Yukawa matrix are predicted to be zero up to corrections from operators of dimension 10 or higher, assuming that the necessary fields required to mediate such off-diagonal operators are even present in the UV theory. The suppressed mixing provides the needed protection against cLFV constraints. A phenomenologically viable realization of the neutrino masses and mixings, on the other hand, requires additional scalars. Thanks to the smallness of the neutrino masses, this can be done consistently without introducing sizeable cLFV.
The 2σ band in parameter space ofL µ−τ that explains (g − 2) µ is denoted with solidblack lines in Fig. 6 (dashed black line denotes the central value), while the colored regions are excluded. The neutrino trident CCFR bound (orange) limits theL µ−τ solution of (g − 2) µ to m X 400 MeV. In Fig. 6 the kinetic mixing parameter was set to ε = g X /10, comparable to the IR contribution from muon and tau running in the loop, cf. Eq. (5.5) below. For larger values of ε, the atomic parity violation constraints (dark purple) become more stringent (see Section 5.1). The upper limit on ε from atomic parity violation makes the CCFR bound very robust, i.e., the trident can not be removed by choosing ε as in the vector category.
The Borexino bound on neutrino-electron scattering (dark green) limits theL µ−τ solution of (g − 2) µ to the X masses above m X 60 MeV. The NSI oscillation bounds, on the other hand, are not relevant for theL µ−τ model, since inL µ−τ the couplings of X µ to electrons are purely axial and, thus, induce only suppressed, spin-dependent effects. Furthermore, inL µ−τ the X boson does not couple to quarks. These two features ofL µ−τ couplings mean that the NSI oscillation bounds are completely avoided in this model.
The constraints from the resonant searches are obtained using the DarkCast code which by default supports only vector couplings. We approximate theL µ−τ bounds with those for a vector model with charges x L 1 ,L 2 ,L 3 = (1, 13 2 , − 13 2 ). For processes where the mass can be neglected, such as NA64, there is no difference between the axial and the vectorial couplings. Also for the BaBar search, the above vector model reinterpretation still approximates rather well the actual bounds [122,123]. Below the di-muon threshold, the dominate decay channel of X is to invisible final states, while above the di-muon threshold there is a sizeable branching ratio for X → µ + µ − . The e + e − decay mode is suppressed because of the charge hierarchy and gives a sub-leading constraint. In all cases, the X decays are prompt in the targeted parameter space. The resonance searches rule out thẽ L µ−τ solution to (g − 2) µ for m X > 2m µ (BaBar) and for m X 100 MeV (NA64).
SinceL µ−τ has chiral couplings, there are also observables that are only important for such cases of combined vector and axial vector couplings to SM fermions, which we discuss next.

Atomic parity violation
Parity-violating atomic transitions can determine whether X couples axially to electrons. For typical momentum transfers in atomic interactions q 2 m 2 X (for the mass range considered here) and the tree-level exchange of X leads to an effective interaction, where the couplings to quarks are due to kinetic mixing. The effect of X on the Hamiltonian of atomic systems is, thus, modeled with a parity-violating point-like interaction between the electric charge of the nucleus and electrons and mimics the parity-violating interaction between the nucleus and electrons mediated by the exchange of the SM Z boson. Hence, the new muonic force can be seen as a modification of the weak charge of the nucleus [124], where Z is the atomic number. The parity-violating transition has been measured to percent-level precision in 133 55 Cs [125], with the latest analysis finding Q exp weak − Q SM weak = 0.65 ± 0.43 [126]. We can use this to place the 95% CL bound on the atomic parity violating (APV) couplings of the X boson, With no particular UV theory in mind, ε is treated as a free parameter, and the APV bound can be evaded by sufficiently reducing the value of ε. 11 Nevertheless, the running of ε can be used to set an approximate lower bound on |ε| in the absence of tuning. Assuming no other BSM particles below the scale of electroweak symmetry-breaking and working in the limit ε g X , the 1-loop running of ε receives its dominant contribution from renormalization group running between the muon and tau mass [62], resulting in For processes with small momentum transfers, we therefore expect ε g X /10. At high scales, where other states, such as the putative S 3 muoquark, are dynamical, these can also contribute to the running of ε, which should be taken into account in any UV completion of the model. The presence of a rather stringent APV bound, Eq. (5.4), has important phenomenological implications. It excludes the possibility of the X coupling to muons being dominated by the kinetic mixing, which could otherwise decorrelate the couplings to muons and muon neutrinos, cf. Section 3.2. This is not possible forL µ−τ . Since in this model the interactions of X µ with nucleus are induced by the kinetic mixing, the COHERENT bound in Fig. 6 (light orange) can be directly compared with the one from APV (dark purple), the latter being stronger.

Electron anomalous magnetic moment
The second observable for which theL µ−τ model is qualitatively different from a vector-like model such as L µ − L τ is the electron anomalous magnetic moment, (g − 2) e . In L µ − L τ the X coupling to electrons is exclusively due to a (presumed) small kinetic mixing with the photon, resulting in a similarly minute modification of (g − 2) e . In theL µ−τ model, on the other hand, the electron carries a nonzero U(1) X charge and receives a much larger contributions to (g − 2) e (the contribution is also numerically larger because of the coupling to electrons is axial, cf. Eq. (2.5)).
The (g − 2) e has been calculated to high precision in the SM [127] but crucially depends on the exact value of the fine-structure constant. Recently, α has been measured precisely in two experiments involving cesium [128] and rubidium [129] atoms, respectively; however, the measurements are internally inconsistent at the level of 5.4σ. With these inputs, the deviation of the measured (g − 2) e [130] from the SM theory prediction is found to be a exp In theL µ−τ model, the axial coupling of X to the electron gives a negative correction to (g − 2) e (see, e.g., Refs. [95,96] A stronger bound is obtained using α from the rubidium measurement, which prefers a positive contribution to (g − 2) e . However, even this constraint is still weaker than the Borexino bound as shown in Fig. 6.

Neutrino-electron scattering
The reactor experiments TEXONO [77] and GEMMA [78] measured the neutrino scattering on electrons and set competitive constraints on the X µ couplings. The relevant bound shown in Fig. 6 is extracted using theν e e − →ν e e − electron recoil energy spectrum from Fig. 16(b) of the TEXONO analysis paper, Ref [77], closely following the procedure described in [131]. ForL µ−τ , we find the 95% CL bounds and 2.8 × 10 −4 < g X × 100 MeV/m X < 3.5 × 10 −4 , where the NP amplitude is roughly twice the negative SM amplitude. The electron recoil energy is 8 MeV, thus, the above EFT limit is valid for the m X range considered in Fig. 6. This bound is only slightly worse than the Borexino bound. The high-energy beam experiment CHARM-II at CERN measured elastic ν µ e − and ν µ e − scattering [79,80]. The target calorimeter was exposed to the horn-focused wide band neutrino beam from the Super Proton Synchrotron with the mean muon neutrino (antineutrino) energy of E ν = 23.7 ± 0.3 GeV (E ν = 19.1 ± 0.2 GeV). The signature of this process is a forward-scattered electron producing an electromagnetic shower measured by the calorimeter. The variable y is defined as the ratio of the electron recoil energy over the incident neutrino energy. Neglecting the electron mass, the differential cross section is given by where where the X charges of theL µ−τ model are given in Eq. (5.1). The unfolded differential distributions in y for ν µ e andν µ e scattering were reported in Fig. 1 of Ref. [79], albeit in arbitrary units. This information was then used in Ref. [79] to determine the SM weak mixing angle by fitting to the shape of the distributions, while requiring no knowledge of the overall normalization. Similarly, inL µ−τ we are able to set an upper bound on g X for given m X by profiling over two arbitrary nuisance parameters describing the absolute normalisations. The excluded region is shown with red color in Fig. 6. This represents the most stringent limit on the model for m X 100 MeV.
The shape analysis is unfortunately ineffective for larger masses where the information on the overall rate is needed. Although the SM prediction was overlaid in Fig. 1 of Ref. [79], the uncertainty on the prediction was not reported. The expected SM cross section depends on the number of factors subject to systematic uncertainties such as the neutrino flux determination. Conservatively assuming the error on the normalisation of the SM prediction to be 30%, i.e., constraining the two nuisance parameters in the fit to be 1.0 ± 0.3, we find the exclusion shown with the red-hatched region in Fig. 6, covering the remaining viable mass window for (g − 2) µ . While the estimate of a 30% error on overall normalization seems plausible to us, if not even conservative (it is an order of magnitude larger than the relative error in the determination of g νe A [80]), it is not sufficient to claim a definitive exclusion of theL µ−τ model. A proper recast of the CHARM-II bound using the detector level events shown in Fig. 1 of Ref. [80] while accounting correctly for the systematic effects could possibly achieve this, but it is beyond the scope of this work.

Possible connections to B-anomalies
We continue by commenting on the possibility that a single vector mediator is behind a simultaneous solution of (g − 2) µ and B anomalies. As shown in Section 4.1, there is a generic upper bound on the vector mass from neutrino trident production and T parameter, m X 4 GeV. In this situation, B → Kνν provides an important bound since it receives contributions from on-shell B → KX decay. Generically, X decays invisibly with a sizeable branching ratio. As shown in Ref. [62], this puts an upper limit on the effective g bs coupling for a given m X . Since the deviation in R K is proportional to the product of g bs and the effective g µµ coupling, the two together imply a lower bound on g µµ for given m X , restricting the viable parameter space to be in the orange region in Figure 1 of Ref. [62]. The range predicted for the effective g µµ coupling from the (g − 2) µ anomaly, on the other hand, does not overlap with the orange band in Figure 1 of Ref. [62], apart from the small region around m X ∼ 2.5 GeV and q A /q V −0.4 (see also, e.g., Ref. [132]). This possibility is ruled out by the constraints on the trident production in the case of a small kinetic mixing. When the requirement of vanishingly small kinetic mixing is relaxed, this conclusion no longer applies in full generality, since the kinetic mixing can be used to remove the trident constraint. One then needs to invoke constraints from resonance searches to rule out explicit models.
A better way to solve the B anomalies in this framework is to extend the field content by a heavy scalar, a TeV-scale muoquark S 3 in the (3, 3, These charge assignments allow a Yukawa coupling (q c S 3 ) of the muoquark to µ but not to e and τ . Furthermore, for x L 2 = −3x q , the dimension-4 diquark operators qqS † 3 are forbidden and proton decay is suppressed. In other words, a non-universal U(1) X enhances the properties of a TeV-scale leptoquark by keeping the accidental symmetries of the SM while still addressing B-anomalies [82][83][84][85][86]. Out of 276 quark-flavor universal U(1) X charge assignments, 273 allow for the above conditions for inclusion of a muoquark [62]. The three models which fail are: the dark photon model, which has x f = 0 for all fermions, and the U(1) B−L and U(1) N i −N j models. The b → sµ + µ − anomaly is resolved by a treelevel exchange of a TeV-scale S 3 muoquark. This exchange leads to an additional V − A contribution to the b → sµ + µ − transitions, as required by data; see, e.g., Refs. [82,83,[133][134][135][136][137][138][139][140][141][142].

Conclusions
A new massive spin-1 boson X µ coupling vectorially to muons, L ⊃ g X q Vμ γ µ µ X µ , can give a one-loop contribution of the right size to explain the 4.2σ discrepancy in (g − 2) µ . Such a bottom-up simplified model is relatively poorly constrained experimentally. The mass of the X boson can be anywhere from 10 MeV (set by cosmological constraints, see Section 2.5) up to 1 TeV (the perturbative unitarity limit, see Section 2.4). Fully exploring this mass window via direct searches may well require a dedicated future collider strategy such as a muon collider [97].
However, additional theoretical inputs more often than not lead to correlations with other phenomenological probes. In general, these either severely constrain the X µ solution to the (g − 2) µ anomaly or predict a signal in the next generation of experiments. For instance, demanding electroweak gauge invariance, the renormalizable couplings of X µ to muons come from two sources: the kinetic mixing, L ⊃ ε Y B µν X µν /2, and the covariant derivative in the kinetic term, which gives couplings of the form L ⊃ g X x L 2L 2 γ µ L 2 + x E 2Ē 2 γ µ E 2 X µ . Each of the two terms comes with a separate set of constraints. Due to electroweak gauge invariance, the covariant derivative term generates couplings of the X boson not just to muons but also to muon neutrinos. This then leads to constraints from neutrino trident production, which is most relevant in the high X mass region. Meanwhile, the kinetic mixing is constrained by electroweak precision data, which is also mostly relevant for heavier X masses starting from around the GeV scale. Even before committing to a specific gauged U(1) X model, the combination of the neutrino trident production and electroweak precision tests, assuming U(1) X breaking by SM gauge singlets, already limits any such solution of (g − 2) µ to relatively light X boson masses, m X 4 GeV, see Section 4.1.
Further phenomenological implications are obtained in specific complete theories that contain the X boson. In this manuscript, we performed a comprehensive survey of spontaneously broken anomaly-free gauged U(1) X models in which the SM matter content is minimally extended by three generations of right-handed neutrinos. Our results are independent of how the U(1) X gauge group is broken, as long as the condensate is neutral under the SM. In the mass window 10 MeV m X 4 GeV we perform a thorough investigation of 419 phenomenologically inequivalent renormalizable models with vector-like charge assignments, allowing for arbitrary kinetic mixing, i.e., for the full set of quark flavor-universal models that have vector-like U(1) X charge assignments for charged SM leptons, and a maximal (finite) U(1) X charge ratio of 10. We find that 7 such models, listed in Table 1, avoid the experimental constraints in a narrow mass window with m X just below 2m µ such that the global tension with all the data, including (g − 2) µ , is less than 2σ. The viable models have charge assignments that are either exactly L µ − L τ or, in most cases, deformations of it in the B − 2L e − L τ direction, which to a large extent avoids the stringent constraints on nonstandard neutrino interactions from neutrino oscillations. The best agreement with the data is obtained for L µ − L τ . The key complementary constraints on the models are due to the searches for nonstandard neutrino interactions, either the bounds from neutrino oscillations or from COHERENT and Borexino.
There are also 21 chiral anomaly-free U(1) X charge assignments with charges in the range [−10, 10]. The contributions to (g − 2) µ are maximized for muon couplings that are as close to vector-like as possible, while the effect of bounds on nonstandard neutrino interactions from neutrino oscillations is minimized for axial electron couplings and for vanishing couplings to quarks. We performed a detailed analysis of a prime candidate of this type, theL µ−τ model, and found that it is excluded in the region relevant for (g − 2) µ . We can reasonably expect that the same is true for the other models with chiral charge assignments.
All of the above U(1) X solutions to the (g − 2) µ anomaly feature lepton flavor-nonuniversal charge assignments. These imply selection rules on the charged lepton mass matrix and force the interaction and mass bases to coincide up to small, potential corrections from higher-dimensional effective operators. This provides a very effective mechanism to suppress charged lepton flavor violation. Indeed, radiative muon and tau decays in general present the biggest challenge for the new physics explanations of the (g − 2) µ anomaly, requiring a rather stringent flavor alignment [44]. In the above models, this protection against such flavor constraints is built into their symmetry structure. Such setups are also interesting in the context of the ongoing B-physics anomalies, allowing for the muoquark mediators to be present [62,[82][83][84][85].

Acknowledgments
We thank Chaja Baruch and Yotam Soreq for many useful discussions. The work of AG and AET has received funding from the Swiss National Science Foundation (SNF) through the Eccellenza Professorial Fellowship "Flavor Physics at the High Energy Frontier" project number 186866. The work of AG is also partially supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, grant agreement 833280 (FLAY). The work of PS is supported by the SNF grant 200020 204075. JZ acknowledges support in part by the DOE grant de-sc0011784.

A NSI oscillation bounds
Here we describe the construction of a χ 2 function that approximates the results of a global fit to the NSI oscillation data [71]. The starting point is an observation that neutrino oscillations in matter depend on the effective parameters 12 (for more details, see e.g. [70,71]) where α = {e, µ, τ } is the neutrino flavor index, Y n ( x) is the neutron-to-proton ratio in matter at position x along the neutrino trajectory, while ε f α are given by 13 where x f,α are the corresponding U(1) X charges, and The U(1) X models we consider have quark flavor-universal couplings, and, thus, the proton and neutron charges are given by the common baryon charge, We make an additional assumption that the approximate χ 2 function can be expressed in terms of the x-independent average E α , which depends on the averaged Y n and is given by We then define the approximate NSI χ 2 function as T denotes the minimum of the χ 2 function.
The components of the covariance matrix C are given by with the standard deviations σ α and the correlation coefficients ρ αβ , which satisfy The parameters η that define the χ 2 function, η = (E 0 e , E 0 µ , E 0 τ , σ e , σ µ , σ τ , ρ eµ , ρ eτ , ρ µτ ) T , (A.9) are then determined from the fit results of Ref. [71], which are provided for various model parameters θ, θ = (x B , x e , x µ , x τ ) T . (A.10) The central quantity that was used to set the bounds in [71] is where χ 2 NSI, SM is the value of the χ 2 function at the SM point. Using the definition in Eq. (A.6), along with ε 0 | SM = 0 and E| SM = 0, we thus have In order to determine the parameters η, we use two quantities that were provided in Ref. [71] for several different models, i.e., for several different values of θ (see the first column in Table 2): • The first quantity is the 2σ bound on g X /m X in a given model, θ i , which is obtained when ∆χ 2 NSI = 4. We convert this to a bound on ε 0 using (A.3), giving theε bnd 0,i values listed in the last column of Table 2. Ideally, these bounds would be reproduced for each θ i by minimizing the ∆χ 2 NSI function, Eq. (A.12), after a judicial choice of the values for nuisance parameters η and Y n . To this end, we solve the equation ∆χ 2 NSI (ε 0 , θ, η, Y n ) = 4 for ε 0 , which results in a function ε bnd 0 ( θ, η, Y n ) , (A. 13) that depends on the model parameters θ, as well as on the χ 2 parameters η, and on Y n . The nuisance parameters need to be chosen so as to minimize the differences betweenε bnd 0,i and ε bnd 0 for each θ i .
• The second set of quantities are ∆χ 2 NSI,min,i , the minimum values of ∆χ 2 NSI for each model θ i , listed in the second column of Table 2. If the nuisance parameters are chosen correctly, ∆χ 2 NSI,min,i should be well reproduced by ∆χ 2 NSI,min ( θ i , η, Y n ) . (A.14) Here ∆χ 2 NSI,min ( θ, η, Y n ) is a function obtained by minimizing Eq. (A.12) with respect to ε 0 .  Table 2. Input data from Ref. [71] used for fitting the parameters of the approximate ∆χ 2 NSI . The value ofε bnd 0,i for θ T i = (1, −2, −1, 0) has been provided by the authors of [71] in private communication in the context of [62].
To find the best values for nuisance parameters we construct a loss function λ, (A.15) where i is an index labelling the different models with parameters θ i for which [71] provides the boundsε bnd 0,i and the minimum values ∆χ 2 NSI,min,i , listed in Table 2. The standard deviations σ ε i and σ χ i are chosen as follows: • In reproducing theε bnd 0,i bounds we allow a nominal 10% uncertainty on the value of the approximate χ 2 function, motivated by how large we expect the deviations between the exact and the approximate χ 2 values to be. We thus set σ ε,i = 0.1ε bound 0,i . (A.16) • The exact χ 2 function used in Ref. [71] to obtain the minima ∆χ 2 NSI,min,i is not Gaussian. We therefore do not expect the approximate χ 2 function to be able to reproduce the values of minima with high accuracy. In fact, we observe that choosing too small values of σ χ,i results in a rather poor fittedε bnd 0,i values. Allowing for very large σ χ,i uncertainties, on the other hand, results in a fit not converging at all. As a compromise, we use the constant values σ χ,i = 2 . (A.17) We minimize the loss function λ using the iminuit [143,144] Python package and obtain a good fit with λ| min = 1. 43 Table 3. Fitted values of the parameters of the approximate ∆χ 2 NSI . the boundsε 0 bound,i with a better accuracy than the assumed tentative uncertainty of 10%. The results of this fit is shown in Table 3. Note that the fitted average neutron-to-proton ratio Y n in Table 3 is very close to the neutron-to-proton ratio averaged over all the Earth Y ⊕ n = 1.051 (but one does not expect to match it exactly, since the oscillation data also include results from oscillation inside the sun).

B Global χ 2 function
In Section 3.3 we use a χ 2 function that combines bounds from various measurements: ∆χ 2 = ∆χ 2 aµ + ∆χ 2 Borex. + ∆χ 2 NSI osc. + ∆χ 2 COHERENT + ∆χ 2 reson. . (B.1) The contribution ∆χ 2 aµ is defined with respect to the ideal NP model, which would be able to exactly reproduce the central value of ∆a µ , Here ∆a µ NP is the shift in the anomalous magnetic moment of the muon due to NP, while ∆a µ exp = 251 × 10 −11 is the difference between the measured value and the consensus SM predictions, with σ(∆a µ ) = 59 × 10 −11 the corresponding error (see also Section [2]). For the SM ∆a µ NP = 0, and thus ∆χ 2 aµ = 17.8. The contribution ∆χ 2 Borex. is due to possible NP contributions to the cross section for 7 Be solar neutrinos scattering on electrons. Averaging the high and low metallicity standard solar model predictions for the solar neutrino fluxes, treating the difference to the two predictions as the systematic error, and adding it in quadrature to the experimental error on the measured Borexino rate, gives R Be = σ Borexino /σ SM = 1.06 ± 0.09 for the ratio of the measured and SM neutrino cross sections. From this we construct where R Be NP = σ SM+NP /σ SM , while R Be exp and σ(R Be ) are the central value and the error of the R Be measurement, respectively. For the SM ∆χ 2 Borex. = 0.44. For NSI oscillations, we use ∆χ 2 NSI described in Appendix A. Since this is only an approximation to the true χ 2 for NSI oscillation bounds, ∆χ 2 NSI can become negative in regions where the approximations in Appendix A become less reliable. We thus take ∆χ 2 NSI osc. = ∆χ 2 NSI if ∆χ 2 NSI > 0, but set ∆χ 2 NSI osc. = 0 if ∆χ 2 NSI becomes negative (signaling the breakdown of our approximations). For the SM ∆χ 2 NSI osc. = 0.
The term ∆χ 2 COHERENT encodes the constraints from the COHERENT measurement. To obtain its value we use the code accompanying Ref. [110], and calculate the χ 2 statistics for the 12T,12E binning. We define ∆χ 2 COHERENT as the difference between the χ 2 value for a particular model and the SM χ 2 value. For SM therefore ∆χ 2 COHERENT = 0. For all the models we require that they pass the bounds on resonance searches as implemented in DarkCast. The term ∆χ 2 reson. is therefore set to zero, if the model passes the NA64 and BaBar 2014 bounds, and it set to a very high value if they do not (in the code we use the numerical value ∆χ 2 reson. = 100 in that case). For the SM ∆χ 2 reson. = 0.