Prospects for a flavour violating $Z^\prime$ explanation of $\Delta a_{\mu,e}$

The apparent tensions emerging from the comparison of experimental data of the anomalous magnetic moments of the muon and electron to the Standard Model predictions ($\Delta a_{\mu,e}$) could be interpreted as a potential signal of New Physics. Models encompassing a light vector boson have been known to offer a satisfactory explanation to $\Delta a_{\mu}$, albeit subject to stringent experimental constraints. Here we explore a minimal extension of the Standard Model via a leptophilic vector boson $Z^\prime$, under the hypothesis of strictly flavour-violating couplings of the latter to leptons. The most constraining observables to this ad-hoc construction emerge from lepton flavour universality violation (in $Z$ and $\tau$ decays) and from rare charged lepton flavour violating transitions. Once these are accommodated, one can saturate the tensions in $\Delta a_{\mu}$, but $\Delta a_{e}$ is predicted to be Standard Model-like. We infer prospects for several observables, including leptonic $Z$ decays and several charged lepton flavour violating processes. We also discuss potential signatures of the considered $Z^\prime$ at a future muon collider, emphasising the role of the $\mu^+\mu^- \to\tau^+\tau^- $ forward-backward asymmetry as a key probe of the model.


Introduction
In recent years, numerous tensions between the Standard Model (SM) expectations and observation have emerged, many (if not most) in relation with high-intensity flavour observables. In addition to the various "anomalous" behaviours associated with B-meson decays, the anomalous magnetic moment of charged leptons (both muons and electrons, albeit the latter to a smaller degree) have been the object of intensive dedicated studies. The anomalous magnetic moment of a charged lepton , defined as allows probing numerous aspects of the SM, and is also instrumental in determining some of its fundamental quantities. Following the disclosed results from the "g-2" E989 experiment at FNAL [1], which are in good agreement with the previous findings of the BNL E821 experiment [2], the current experimental average for the muon anomalous magnetic moment [1] is given by which should be compared to its SM expectation. Prior to the most recent lattice QCD based computation of the hadronic vacuum polarisation contribution by the BMW collaboration [3], the SM prediction -as compiled by the "Muon g − 2 Theory Inititative" [5] (see also [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] -was found to be a SM µ = 116 591 810 (43) × 10 −11 .
When compared, the experimental average and the latter SM prediction lead to the following 4.2 σ tension between theory and observation 1 Under the assumption of a significant tension between theory and observation, as given by Eq. (4), the need for new physics (NP) capable of accounting for such a sizeable discrepancy is manifest; several minimal, as well as more complete NP models, have been thoroughly explored in the light of recent experimental results (for a review see, for example, [26] and references therein). The magnetic moment of the electron has also been at the origin of possible new tensions, upon comparison of the experimental value to the SM prediction (depending on the value of α e that is used for the computation of the latter): using α e as extracted from measurements using Cs atoms [27,28], one is led to corresponding to a −2.5 σ deviation; in [29], a more recent estimation of α e was obtained, this time relying on Rubidium atoms, and the new determination of α e (implying an overall deviation above the 5 σ level for α e ) now suggests milder discrepancies between observation and theory prediction, ∆a Rb e = 0.48 (0.30) × 10 −12 , corresponding to O(1.7 σ) deviation. Other than (possibly) signalling deviations from the SM expectation, it is interesting to notice the potential impact of both ∆a Cs e and ∆a µ : other than having an opposite sign, the ratio ∆a µ /∆a e does not exhibit the naïve scaling ∼ m 2 µ /m 2 e (expected from the magnetic dipole operator, in which a mass insertion of the SM lepton is responsible for the required chirality flip [30]). This behaviour renders a common explanation of both tensions quite challenging, calling upon a departure from a minimal flavour violation (MFV) hypothesis, or from single new particle extensions of the SM (coupling to charged leptons [31][32][33]). Notice that the pattern in both ∆a Cs e and ∆a µ can be also suggestive of a violation of lepton flavour universality (LFU).
Several new physics constructions have been put forward to simultaneously explain the tensions in a e and a µ (see, for example, ). Among the most minimal models, extensions of the SM gauge group via additional U (1) groups have been intensively explored, as these offer several appealing features. The new (vector) boson (as well as other potentially present state) and new associated neutral currents can open the door to extensive implications both for particle and astroparticle physics, across vast energy scales [60]. In addition to their potential in what concerns well-motivated dark matter mediators see, for instance [61][62][63][64], Z extensions of the SM have been considered in the context of flavour physics, in particular in what concerns B-meson decay observables (especially b → s transitions) [61,[65][66][67][68][69][70][71][72][73][74][75][76][77][78][79][80].
The most minimal SM extensions leading to a (light) Z boson have mostly been oriented towards constructions featuring flavour-conserving couplings of the new vector boson to matter (quarks and leptons). Models in which the Z couples to hadrons have a severely constrained parameter space due to conflict with numerous bounds from meson decays and oscillations, as well as from atomic parity violation experiments [81][82][83]. Likewise, and despite their potential to address the B-meson decay "anomalies" (and/or tensions in a µ ), lepton flavour conserving (although not necessarily lepton flavour universal) leptophilic Z extensions have also been strongly constrained, as a consequence of the extensive implications for (lepton) flavour observables, high-and low-energy neutrino observables, LHC phenomenology, as well as dark matter direct detection experiments [65,77,.
A NP construction relying on the addition of a new vector boson usually calls for the extension of the SM gauge group by at least an additional U (1). Likewise further states must be added, for instance to explain the origin of a massive vector or to provide a rationale for the peculiar flavour structure. However, as a first step to evaluate the phenomenological viability of SM extensions via leptophilic flavour violating (light) Z bosons, one can rely in a minimal, ad-hoc approach. In this work, we consider a toy-model construction in which only new interactions with the leptons are added to the SM Lagrangian, allowing for left-handed and right-handed couplings (g αβ L,R¯ α β Z ). No assumption is made on the underlying gauge group, nor on the actual mechanism leading to its breaking. Motivated by an explanation to ∆a e,µ , we only hypothesise the existence of a massive (albeit light) vector boson, and investigate the constraints on its g αβ L,R couplings: in addition to addressing whether or not it can simultaneously saturate the tensions for both the electron and the muon anomalous magnetic moments, we consider numerous constraints arising from electroweak precision observables (EWPO), rare (lepton) flavour violating transitions and decays, as well as LFU tests in Z and tau-lepton decays, providing in most cases a full computation (beyond renormalisation group (RG) approximation) of the relevant quantities (see also [106][107][108][109] and [110], albeit in the context of effective theory studies). As we proceed to discuss, the latter bounds lead to very stringent constraints on the Z parameter space: perturbativity of the couplings allows inferring an upper limit for its mass (m Z ), and the flavour constraints on the couplings suggest that ∆a e (and ∆a τ ) should be SM-like. Interestingly, several charged lepton flavour violating (cLFV) observables could be within future experimental sensitivity.
Our analysis is further complemented by investigating the prospects of such a leptophilic Z in what concerns a future muon collider. Muon colliders have received increasing attention in recent years [111][112][113][114] and offer promising testing grounds for beyond the SM (BSM) constructions with preferred couplings to leptons. In line with recent studies (see, e.g. [91]), in this work we also discuss the t-channel Z contribution to di-tau pair production, µ + µ − → τ + τ − , comparing the deviations with respect to the SM expectation regarding the production cross section, and the associated forwardbackward asymmetry (A FB ). As discussed here, our findings suggest that a muon collider could indeed offer good testing (and discovery) grounds for a leptophilic Z capable of explaining the tension in ∆a µ , or conversely being capable of falsifying the proposed explanation for the tension in the anomalous magnetic moment of the muon.
We emphasise that in this work we have focused on studying the low-energy effects induced by the presence of a new neutral vector boson exhibiting distinctive couplings exclusively to the lepton sector. This is but a first step towards the model-building of ultraviolet (UV)-complete SM extensions whose particle content would feature such a neutral vector boson, in addition to other new states (scalars, fermions, ...).
The manuscript is organised as follows: in Section 2, we describe the minimal (ad-hoc) construction we will explore, subsequently describing the most stringent constraints on flavour violating Z couplings to leptons in Section 3, stemming from complying with relevant bounds -from both cLFV transitions and from observables sensitive to lepton flavour universality violation (LFUV). In Section 4 we address the possibility of saturating ∆a e,µ while complying with the considered bounds. In Section 5 we discuss the impact that the latter constraints might have on the prospects for the maximal predictions of a number of charged lepton flavour violating observables. This is followed by a discussion of the prospects of such a state at a future muon collider (Section 6). We then present our summarising conclusions and a brief overview.
2 Leptophilic cLFV Z : a minimal model As mentioned in the Introduction, we consider here a minimal extension of the SM via a single neutral vector boson Z (without specifying the underlying gauge group). Generic Z extensions of the SM, especially those in which the new mediator is lighter than the electroweak (EW) scale, Λ EW , are subject to stringent constraints, both from direct searches and from electroweak precision observables.
Beyond the SM constructions in which the Z couples to hadrons are constrained from light meson decays, π 0 → γZ (Z → ee) and K + → π + Z (Z → ee), as searched for at the NA48/2 [81] experiment, as well as by searches for φ + → η + Z (Z → ee) at KLOE-2 [82]. Likewise, rare meson decays (e.g. B (s) → + − ), neutral meson oscillations (K 0 −K 0 and B (s) −B (s) mixing), and atomic parity violation [83], all play an important role in constraining hadrophilic Z extensions of the SM. Most of these constraints can be avoided by considering specific extensions in which the Z only couples to the lepton sector (neutrinos and charged leptons) 2 . Nevertheless, flavour conserving Z couplings to leptons, i.e. Z α α , also lead to NP contributions that are potentially in conflict with numerous observables: notice that the non-observation of a Z at electron beam dump experiments (SLAC E141, Orsay, NA64 [115,116]), in dark photon production searches (KLOE-2 experiment [82], BaBar [117]) or at parity-violation experiments (SLAC E158 [118]), set severe bounds on Z ee couplings. The other diagonal couplings (second and third generation, as well as to left-handed neutrinos) are also constrained byν e -electron and ν µ -electron scattering, in particular from the data of the TEXONO [119] and CHARM-II [120] experiments, respectively. Diagonal couplings to muons are further severely constrained from neutrino trident production (ν µ N → ν µ N µ + µ − ) [121].
In order to evade these severe bounds one can envisage NP models in which the new Z has no flavour conserving couplings to leptons; in our analysis, we thus consider a strictly leptophilic flavour violating (light) Z , leading to the presence of the following new terms in the interaction Lagrangian: In the above, the indices run over α, β = e, µ, τ , with α = β; P X = P L,R are the chiral projectors and g αβ X = g αβ L,R are the new couplings hermitian matrices. Notice that throughout the study, and for simplicity, we will always consider real couplings. In the "toy-model" here considered, we do not extend the lepton sector to account for massive neutrinos; only left-handed (LH) neutrinos are present and, due to SU (2) L gauge invariance, they couple to the Z through the same LH charged lepton couplings, leading to 2 Notice that higher order processes -typically at the one-loop level -, arising from gauge kinetic mixing can still lead to new contributions to the "hadronic" observables. Furthermore, as we proceed to discuss, we will only consider flavourviolating couplings to leptons and thus gauge-kinetic mixing between the photon and the Z only arises at the two-loop level. Kinetic mixing ( kin ) is therefore strongly suppressed, of the order kin eg eµ g eτ g µτ /(256π 4 ) log(µ 2 /m 2 Z ). In what follows, these effects will not be taken into account.
Despite the simplicity of this BSM realisation, it is clear from the structure of the Lagrangian in Eq. (8) that the leptonic couplings of the new neutral gauge boson can potentially lead to extensive contributions to several leptonic and EW precision observables -including the desired sizeable  contributions to the charged lepton anomalous magnetic moments. Before proceeding to discuss this in the following sections, a few comments are in order concerning formal aspects related to possible UV completions for this ad-hoc phenomenological construction, especially regarding anomaly cancellations. Without attempting a full discussion, let us briefly consider the new gauge anomalies emerging in this context, possibly in association with the presence of a new U (1) . First, the simplest ones are the Z − W 1,2 − W 1,2 and Z − B − B one-loop triangle diagrams involving SU (2) L and U (1) Y gauge fields. The vanishing of these anomalies respectively requires Tr(g L ) = 0 and Tr(g R ) = 0. Our peculiar hypothesis of purely off-diagonal leptonic Z couplings provides the simplest solution, which in turn also partly justifies this assumption. Concerning the vanishing of the Z − Z − Z triangle anomaly, this requires 2Tr(g 3 L ) = Tr(g 3 R ). This can be achieved by setting one of the three off-diagonal couplings of each chirality to zero, as this trivially ensures the equality of the traces (identical to zero). This choice can also further prevent the appearance of flavour-conserving Z vertices at the one-loop level. Finally, the Z − Z − W 3 triangle is proportional to Tr(g 2 L ) (a positive definite quantity), and also to the vanishing sum of the lepton and neutrino isospins. The only difficulty then lies in association with the triangle Z − Z − B anomaly (or equivalently Z − Z − γ), whose vanishing would require Tr(g 2 L ) = Tr(g 2 R ), in conflict with the results of the subsequent phenomenological analysis 3 . New fields should thus be added: on top of a new scalar responsible for symmetry breaking and for giving a mass to the Z boson, additional (heavy) charged fermions will be needed to properly cancel the above gauge anomaly, and to provide a realistic origin of the (strictly) flavour violating couplings.

Constraining flavour violating Z couplings to leptons
In addition to the stringent constraints that have been obtained for the lepton flavour-conserving (LFC) couplings of a new Z boson [82,[115][116][117][118][119][120][121], there is an extensive array of bounds on its flavourviolating couplings to leptons, as a consequence of the associated experimental limits. In particular, the bounds on the g αβ L,R couplings (with α = β) stem from both lepton flavour violating transitions and decays, as well as processes sensitive to the breaking of LFU (the latter being an indirect consequence of the former in the present case). Thus, and before addressing the ultimate requirement that should be fulfilled by this class of simplified leptophilic Z extensions (i.e. the tensions in (g − 2) ), in this section we address the constraints on g αβ L,R arising from cLFV and LFUV limits, including Z decays, as well as several leptonic processes.

Constraints from Z decays
The presence of the new Z can lead to new (higher order) contributions to both lepton flavour conserving and lepton flavour violating (LFV) Z decays 4 .
Concerning flavour conserving Z decays, the new Z -mediated loops (and the non-negligible interferences between the SM-like processes and the NP ones) will lead to new contributions to the individual decay widths and, as expected, to modifications of the effective Z couplings. In order to assess the impact of the current experimental measurements, we consider the following ratios of decay widths, which further have the advantage of allowing the cancellation of quantum electrodynamics (QED) corrections in the theoretical predictions. Moreover, these ratios are also probes of the lepton flavour universality of Z-boson couplings. We thus consider the ratios , with α = β = e, µ, τ .
The SM prediction for these ratios (at 2-loop accuracy) are [125] with negligible associated uncertainties. These should be compared with the corresponding experimental values [126], As can be seen, the experimental measurements are in good agreement with their respective SM predictions, thus placing strong bounds on any NP contribution. In what concerns the NP contributions, we estimate the modified individual partial widths as where Γ SM full is given in [125] at 2-loop accuracy, Γ SMtree−Z is the interference term between the SM tree-level contribution and the Z -mediated 1-loop diagrams and Γ Z purely consists of the Z -mediated 1-loop diagrams. In particular, the interference term between the SM tree-level diagrams and the Z contributions is at the source of important corrections to the Z-boson leptonic partial widths, which in turn leads to stringent constraints on the Z couplings to SM leptons. The amplitude of the Z contributions is given by where u α , v α denote the α spinors, * µ (q) the Z-boson polarisation vector, q is the Z momentum and C V,T X are the vector and tensor coefficients computed from the decay "triangle" diagram (i.e. oneloop vertex correction), which are given in Appendix A.1. (These calculations were done relying on Package-X [127].) In view of the structure of the couplings of the Z to leptons -cf. Eq. (8) -, even in the absence of a specific mechanism of neutrino mass generation, new contributions to the so-called invisible Z decay width are expected: these are fuelled by the non-vanishing Z ν α ν β couplings (∝ g αβ L ). However, due to the observed smallness of m ν i (or zero, in the strict SM limit of massless neutrinos) the Z -mediated 1-loop contributions turn out to be always negligible. Moreover, notice that the Z can also mediate loop decays leading to cLFV Z-decays 5 , Z → ± α ∓ β . The relevant expressions of the associated cLFV Z decay amplitude are collected in Appendix A.2. The current limits from ATLAS on Z → e ± µ ∓ [128], OPAL on Z → e ± τ ∓ and on Z → µ ± τ ∓ [129] (see Table 1) are not sufficiently strong to constrain the different Z couplings; in fact, it turns out that the NP contributions to the decay rate are very small: at most one has BR(Z → ± α ∓ β ) ∼ O(10 −14 ) (throughout the model's parameter space which will be numerically explored in subsequent sections).

Leptonic processes: rare decays and transitions
Lepton decays offer numerous probes of the potential contribution of a light leptophilic Z , and lead to extensive constraints on its flavour-violating couplings. Among these processes one finds cLFV muon and tau three-body decays (µ → 3e, τ → 3e, τ → 3µ, τ − → µ − e + e − , τ − → µ − e + µ − , τ − → e − µ + µ − and τ − → e − µ + e − ), most of them mediated by Z and photon-penguin, and radiative decays (µ → eγ, τ → eγ and τ → µγ). As mentioned in the Introduction, we present here a full calculation of the one-loop contributions to LFV operators beyond leading-log approximations (see e.g. [109]) 6 . In the muon sector, the Z can also mediate Muonium oscillations as well as µ − e conversion. Finally, one must also consider the impact of the Z on τ → νν decays, and on the associated ratios of decay rates.

cLFV 3 body decays: α → β γ δ
Depending on the mass of the decaying lepton, and as illustrated in Fig. 1, several final-state flavour configurations are possible. The (effective) Lagrangian governing these transitions can be written as [146] in which L γ includes photon exchanges (dipole and anapole), and L 4 encodes the interactions relevant for Z-penguins and tree-level decays. The former (L Z ) can be cast as with {X, Y } = {L, R} and Figure 1: Feynman diagrams contributing to cLFV three-body decays: from left to right, photon and Z penguins, and tree level Z exchange.
In the above, the vector and tensor form factors, F V,T X , are computed from cLFV Z decays in the limit of vanishing external fermion masses (m α , m β → 0) and of vanishing momentum transfer (i.e. q.q → 0); moreover one hasX = R, L (for X = L, R). One thus finds is minimally substracted (MS-scheme) and the remnant 't Hooft scale set to µ 2 = m 2 α . The tree-level interaction (encoded in L tree ) is given by One thus finds for the "effective" 4-lepton interactions: in the above notice that F Z,αβ XY only contributes in the case of γ = δ see Eq. 15. Concerning the photon-exchange contributions to the cLFV 3-body decays, one has in which e is the electric charge, q denotes the photon 4-momentum and the anapole coupling In turn, this leads to the following photon-penguin amplitude: with Q δ the electric charge of δ . Taking for simplicity (and clarity) the limit of vanishing external lepton masses, m α , m β → 0 (notice however that we keep all terms in the numerical computation of the decay rates), the anapole form factors entering the above equation are given by in which i denotes the internal lepton, x i = m 2 i /m 2 Z , and with f (x) defined as In the limit m β , q 2 → 0 (again, we emphasise that in the numerical computations we do take into account m β = 0), the dipole coefficients are given by where (23)). After considering the most general case for the 3-body decays, we proceed to address some simpler realisations, α → β β β and α → β γ γ , α → γ β γ .

cLFV 3-body decays: α → β β β
This subset of processes receives contributions from Z-and photon-penguins, as can be seen in Fig. 1(a). The associated branching ratio is given by: (27) in which τ denotes the lifetime of the decaying fermion, and where we have introduced the vector contribution, The relevant form factors have been defined in the previous subsection.

cLFV 3-body decays: α → β γ γ
Another relevant case is that of α → β γ γ (with β = γ), corresponding to the processes τ → eμµ and τ → µēe. In the model under consideration these decays can occur at loop-level (penguin diagrams) and at tree-level, see Fig. 1(b). The corresponding branching ratios are given by: As before we define distinct vector contributions (for XX = LL, RR): that include tree-level, Z-and γ-penguin mediated exchanges, and which does not include any tree-level contributions 8 .

cLFV 3-body decays: α → γ β γ
Finally we consider the decay α → γ β γ (with β = γ), that can only arise from tree-level Z exchanges, see Fig. 1(c). The physical processes corresponding to this type of decay are τ → eμe and τ → µēµ, and the corresponding branching ratio is given by: 3.2.5 cLFV radiative decays: α → β γ These correspond to higher order processes, mediated by the light Z in the loop, and the associated branching ratios are given by [149] BR in which all the quantities have been already defined.

Neutrinoless muon-electron conversion in nuclei
The cLFV Z interactions can also lead to contributions to µ − e transitions in muonic atoms. Following [150], the (coherent) conversion rate can be cast as with D, V (p) , V (n) the overlap integrals (see [150]). Here Γ capt is the muon capture rate, and G F is the Fermi constant. The relevant quantities entering the above equation are given by with K µe 2L/R defined in Eq. (26) and In the above, Q q is the electric charge of the quark q, and with F V,T V X given in Eq. (18) and we clarify that in the above g Z,q V denotes the vector coupling of the Z-boson to a quark q. Figure 2: Feynman diagrams for Muonium-antimuonium oscillations (tree-level Z mediated, s and t channel exchanges).

Muonium oscillations
Another relevant observable to consider are spontaneous Muonium-antimuonium oscillations. Muonium is a bound state of an electron and an antimuon, Mu = µ + e − . In presence of cLFV interactions, this system can oscillate into antimuonium Mu = µ − e + . In the model under consideration Mu − Mu conversion can occur at tree-level (as shown in Fig. 2), from both s− and t-channels Z exchange. The relevant effective Lagrangian for the oscillations can be cast as [151]: where Q i are the four-fermion operators responsible for the Mu − Mu transitions Moreover, notice that the presence of an external magnetic field also induces an energy splitting, further contributing to Mu − Mu mixing. In the present LFV leptophilic Z model, the transition probability can be written as where |c F,m | 2 denote the population of Muonium states and the factor X encodes the magnetic flux density [151]. Notice that in the above we have only considered vector couplings (so that Q 4 , Q 5 = 0). For Hermitian couplings g eµ L,R (we recall that in our study we work with real couplings), one finds The new contributions to the transition probability P (see Eq. (39)) should be compared with the bound set by the PSI experiment [152]

Constraints from LFU in τ decays
The comparison of the SM-allowed τ → ν τ ν decays (mediated at tree-level by W -boson exchange) allows to construct the following ratio, which can be used as a probe of LFU in lepton decays: 8 The latter will be taken into account in the scalar contributions of the type A S XY = −2F αγγβ XY that stem from the Fierz-transformed tree-level operators to match the penguin ones.
The measured τ leptonic decay ratio R τ µe | exp from ARGUS [153], CLEO [154] and BaBar [155] is found to be consistent with the SM prediction R τ µe | SM = 0.972564 ± 0.00001 [156]; the HFLAV collaboration reports the following global fit [157]: Experimentally, it is not possible to disentangle the different neutrino flavours, so any LFV tree-level decay (τ → ν γ ν δ , with γ, δ = e, µ, τ ) will contribute to this observable; tree-level contributions from the new light Z can thus potentially compete with the SM processes, leading to deviations from the observed R τ µe | exp . The relevant contact interactions at the source of the new contributions to the decays (generically cast for α → β ν γνδ ) can be written as: where we have introduced C after Fierz-transforming the SM charged current operator, one has The SM contribution, only defined for the process α → β ν ανβ , reads with G eff F the effective Fermi constant in the presence of NP in µ → eν e ν µ decays; in fact, and as a consequence of the modification of the muon lifetime, a comparison of the new rate of µ → eν e ν µ decays with the SM prediction allows defining G eff F as .
The tree-level NP contribution to the decay α → β ν γνδ is given by thus leading to the following expression for the τ → νν decay width and the associated functions defined as Further QED corrections and effects from the non-local structure of the W -boson propagator are included in r αβ RC 9 : where α e (m α ) is the running electromagnetic fine-structure "constant" (at m α ).
Generally, the above corrections to the Fermi constant G F due to the presence of New Physics in the Michel decay of the muon should be propagated to other electroweak observables such as the invisible Z-decay and the mass of the W -boson. However, as discussed in detail in the following sections, experimental data constrains the g eµ couplings to be very small and thus corrections to G F are negligible for all practical purposes.

Implications for the anomalous magnetic moments
We now address whether or not the model under consideration (a light, leptophilic, cLFV-interacting Z boson) allows to account for the observed a e,µ , while complying with the numerous bounds from boson decays as well as rare leptonic processes presented in the previous section. We first discuss the new contributions mediated by the (light) Z , and then (numerically) explore the model's parameter space (spanned by new boson's mass, and by its left-handed and right-handed cLFV couplings, g αβ L,R ). The new physics contribution to a from the additional Z boson arises at one-loop (see Fig. 3) and, following [159][160][161], it can be expressed as: with g V , g A = (g L ± g R )/2 denoting the vector and axial-vector couplings. The sum runs over the internal leptons (whose flavour is necessarily distinct from that of the external legs), and the loop function F (λ, i ) is defined as: in which i = m i /m , m i being the mass of the internal fermion and λ = m /m Z . Notice that the opposite sign between vector and axial-vector loop functions leads to partial cancellations, thus potentially opening the door for an explanation of both ∆a µ and ∆a e . Throughout this section, in what concerns constraining the model's parameter space, we will fix the Z mass to m Z = 10 GeV (unless otherwise stated); this is a simplifying approach, which allows to evade the numerous constraints which would otherwise arise from the extensive searches at low-energy experiments such as indirect signals at B-factories [162], and flavour violating τ → µZ decays [106][107][108]. At the end of the discussion (Section 4.4), we will revisit this working hypothesis, and explore a wider range for m Z .
Before carrying out the numerical study, let us briefly clarify some points regarding the associated statistical analysis. Leading to the (best-fit) contours presented in several of the plots of this section, we note here that contours of a single observable correspond to model predictions within the 1 σ region of the experimental measurements; for combinations of observables (e.g. the Z-decay LFU ratios and the global contours), we construct a combined likelihood as a product of the experimental probability distribution functions (pdf) of the observables O i , evaluated at a point of the model's predictions given a set of input parameters p: For the observables here considered, theoretical uncertainties due to SM input parameters are negligible and we approximate the observables' probability density functions as uncorrelated gaussians. Thus, the relative (log)-likelihood function can be approximated as a simple χ 2 -function given by in which O exp i and σ exp i respectively denote the mean value and the 1 σ (gaussian) uncertainty of the O i observable measurement. The χ 2 -function is then minimised in terms of the model's parameters.
We further consider the test statistic ∆χ 2 ≡ χ 2 ( p) − χ 2 min in order to establish confidence intervals; in particular we assume that this quantity is distributed as a χ 2 random variable with n = 2 degrees of freedom. The combined contours then correspond to the k σ thresholds where the 2-dimensional cumulative χ 2 distribution reaches the probability P k σ (for 1 σ ∆χ 2 ≈ 2.3, for 2 σ ∆χ 2 ≈ 6.2), defined as the probability for a gaussian random variable to be measured within k standard deviations from the mean.

Accommodating the anomalous magnetic moment of the muon
We begin by discussing the regimes for the Z α β couplings (i.e. g αβ L,R ) which allow to alleviate the current tension in ∆a µ , see Eq. (4), taking into account the impact regarding the numerous cLFV and LFU observables discussed in Section 3. The NP contributions to (g − 2) µ can be associated with the exchange of either an electron or a tau lepton in the loop. The first possibility is strongly disfavoured, as it would lead to a negative shift in ∆a µ . Let us then consider the second possibility: in Fig. 4, we present two representations of the plane spanned by the µ − τ couplings, (g µτ L − g µτ R ). As already mentioned, we set m Z = 10 GeV. On the left panel of Fig. 4 we display the ∆a µ -favoured regimes relying on µ − τ cLFV Z couplings (all others being set to zero), as well as the most important associated constraints, which in this case arise from conflicts with the LFU-probing ratios, R τ µe and R Z αβ . Notice that the constraining role of R τ µe is much more important than that of R Z αβ ; the latter can be easily accommodated in wide regions of the parameter space (the only exception being large values of g µτ L,R 3 × 10 −1 ). The other cLFV observables discussed in Section 3 play a far less restrictive role, and we do not discuss them here. For completeness, let us notice that corrections to the invisible Z decay width due to one-loop contributions are negligible, ∆Γ(Z → inv.) 1 keV, throughout the model's parameter space.
As can be seen, all experimental observations (i.e. ∆a µ and the LFU constraints) can be accommodated for g µτ L ∼ O(10 −3 ), and g µτ R ∼ O(10 −2 ). For clarity, this region has been expanded in the right panel of Fig. 4, where we have also displayed the corresponding regimes with negative values of the couplings (notice however that one must have same-sign couplings). In addition to ∆a µ and R τ µe , we also display the 1 σ and 2 σ contours of the global fit to the observables. Finally, also visible from the right panel of Fig. 4 is the apparent relation emerging for the left-and right-handed Z couplings: the best-fit regions to saturate ∆a µ lie around g µτ R ∼ O(15) × g µτ L .

Accounting for the anomalous magnetic moment of the electron
We now consider how the new Z contribution(s) can further help addressing the recently identified tension in ∆a Cs e (and ∆a Rb e ), firstly studying it independently of ∆a µ . In Fig. 5 we thus display the viable regimes in the plane spanned by e − µ couplings, (g eµ L − g eµ R ), again setting m Z = 10 GeV (and fixing the remaining couplings to zero). As can be readily seen from the left panel, the current tension arising from the determination of α Cs e cannot be accounted for with same-sign couplings due to conflicts with the experimental bounds on LFU Z decays (i.e. R Z αβ ). While the latter bounds still allow for a light Z explanation of ∆a Rb e , the preferred regime of g eµ L,R is incompatible with the current upper bounds on cLFV Muonium oscillations. The most promising regions identified in the left panel of Fig. 5 (logarithmic scale) are presented in detail in a linear scale in its right panel, in which we omitted the R Z αβ constraints (but highlight that these are indeed satisfied in all the regions displayed). Although we will return to this in the following subsection, the regimes of g eµ L,R favoured by an explanation of ∆a e would also contribute to the muon (g − 2), in fact leading to a negative shift in its value, and thus worsening the discrepancy with respect to the SM prediction.  Figure 4: Constraints on Z couplings, g µτ L,R , and prospects for (g − 2) µ . On the left, allowed regimes for ∆a µ in the same-sign (g µτ L -g µτ R ) parameter space, in agreement with the experimental bounds on LFU Z decay ratios (light orange) and R τ µe (purple); the dark orange region denotes the 1 σ ∆a µ favoured regime. On the right, detailed (linear) view of the viable g µτ L vs. g µτ R parameter space, with the dashed curves now corresponding to the 1 σ and 2 σ global fits (see text).  Similar to what was done for (g − 2) µ , one can also rely on Z − eτ couplings (i.e. g eτ L,R ); the prospects for ∆a Cs, Rb e are shown in Fig. 6, which displays a view of the (g eτ L − g eτ R ) plane, with all other couplings set to zero. As clearly seen from the plots, in order to saturate ∆a Cs e the couplings must have opposite signs (or one of them be zero); analogous to what was observed for ∆a µ , ∆a Rb e can be explained relying on a combination of same-sign g eτ L,R couplings. Although all relevant cLFV and LFU observables were taken into account, the most constraining limits arise from R τ µe . Interestingly, notice that such non-vanishing values of g eτ L,R also lead to contributions to ∆a τ ; these are negative, and typically O(10 −7 ) (in fact much smaller than the uncertainty associated with the SM prediction for (g − 2) τ ).

A joint explanation for ∆a µ and ∆a e ?
After having independently considered the new contributions to the muon and to the electron anomalous magnetic moments (taking non-vanishing individual couplings at a time), one must now address the possibility of a joint explanation to both tensions. In order to do so, and in view of the very restricted regimes for the g µτ L,R couplings which allowed explaining ∆a µ , we thus fix g µτ L,R as to comply with (g − 2) µ (we set g µτ L = 0.0024 and g µτ R = 0.036, see Fig. 4), and vary g eµ(τ ) L,R . The prospects for a joint explanation of ∆a µ and ∆a e are depicted in Fig. 7: on the left (right) panel we explore the viable regimes for the g eµ L,R (g eτ L,R ) couplings. In addition to the constraints already identified in Figs. 5 and 6, the simultaneously presence of g eµ L,R = 0 (or g eτ L,R = 0) and g µτ L,R = 0 opens the door to sizeable contributions to cLFV processes: in the left panel of Fig. 7, it is visible how the previously most stringent process (i.e. Muonium oscillations) is now strikingly superseded by various cLFV tau decays. In particular, τ → eμµ and τ → µēµ directly exclude any regime with g eµ L,R 10 −5 . Similar results are obtained when jointly considering the effects of g eτ L,R = 0 and g µτ L,R = 0: as displayed in the right panel of Fig. 7, the constraints arising from rare cLFV muon decays are significantly more restrictive than those arising from R τ µe . Three-body muon decays already constrain g eτ L(R) 10 −6(−5) , with the muon cLFV radiative decay further imposing g eτ L(R) 10 −8(−7) . In summary, for m Z = 10 GeV, an excellent fit to the data (saturating ∆a µ , and complying with R Z αβ , R τ µe ) can be found for the bounds from τ → µēµ and µ → eγ then respectively imply the following upper limits for the combinations of couplings, The above discussion strongly suggests that a joint explanation 10 to the tensions in the light charged leptons anomalous magnetic moments is clearly precluded in view of the cLFV constraints.
Should the current hypothetical model of a light Z with flavour violating couplings to charged leptons be indeed an explanation to ∆a µ , then one is led to mostly SM-like scenarios to both (g − 2) e and (g − 2) τ . On the left, allowed regions in the (g eµ L -g eµ R ) parameter space, in agreement with constraints from cLFV τ decays (τ → 3e in brown, τ → eγ in orange, τ → eμµ in light blue and τ → µēµ in pink) and Muonium oscillations (blue); the red and the green regions respectively denote the 1 σ (g − 2) Cs e and (g − 2) Rb e favoured regimes. On the right, allowed regions in the (g eτ L -g eτ R ) parameter space, in agreement with R µe (τ → νν) (purple), µ → 3e (brown), µ → eγ (orange) and µ − e conversion in Gold (turquoise); the red and the green regions respectively correspond to the 1 σ (g − 2) Cs e and (g − 2) Rb e favoured regimes.

Heavier mediator regimes
For simplicity, we have so far considered an illustrative value for the mass of the Z . In order to complete the study, it is important to discuss to which extent the conclusions of the previous subsections hold for other regimes of m Z . We thus explore regimes for the Z couplings allowing to account for ∆a µ (i.e. g µτ L,R ), for varying masses of the mediator. Relying on the findings of Section 4.1, we set g µτ R ∼ 15 × g µτ L , now for m Z ∈ [2 GeV, 1 TeV] (the lower bound allowing to escape the constraints that emerge for very light states below the m τ threshold, especially τ → µZ [106][107][108], as discussed at the beginning of this section). In Fig. 8 we display information similar to that summarised in Fig. 4, presenting the best global fit regions (1 and 2 σ) in the model's parameter space now spanned by g µτ R and m Z . As before, we colour-code the regimes in agreement with the most constraining bounds (for g µτ L,R couplings, the LFU ratios R τ µe and R Z αβ ). The results displayed suggest that light Z mediators, with masses m Z ∈ [10 GeV, 200 GeV], with associated couplings 0.01 g µτ R 1 offer the best prospects to account for the current tensions in ∆a µ .
Another insight on the best-fit regimes of the model's parameter space, favoured by an explanation to (g − 2) µ , is presented in Fig. 9, in which g µτ L,R are independently fitted for varying values of m Z (for details on the fit procedure, see the beginning of this section).
On the left panel of Fig. 9, we present the results of the best-fit for g µτ L,R taking m Z to lie in the range [2 GeV, 200 GeV]. As can be seen, and despite having taken independent input values for the couplings, the fit leads to a correlation between g µτ L and g µτ R , strengthening the original findings (and underlying assumption leading to Fig. 8). For completeness, we present in the right panel of Fig. 9 a similar study, but now including a much wider interval for m Z , up to 10 TeV. However, for mediator masses O(1 TeV), one clearly runs into a scenario for which no satisfactory fit can be achieved (as the couplings become non-perturbative); in turn this allows to infer a "soft" limit for the validity of the model, m Z ∼ O(1 TeV). g µτ L ± 1 σ g µτ R ± 1 σ g µτ L,R > √ 4π Figure 9: Fits of g µτ L,R as a function of m Z . In blue (orange), 1 σ estimate of g µτ L(R) ; the dashed lines denote the best fit value of the couplings for a given m Z . On the right panel, regimes for which the fit is not satisfactory (non-perturbative couplings) are represented in grey.

Prospects for cLFV probes
Interestingly, one can also estimate the (theoretical) upper bounds for several observables, arising from regimes of g µτ L,R and m Z saturating ∆a µ , while in agreement with all other constraints. The most relevant predictions are collected in Fig. 10, in which we display the projections for the maximal values of several observables as a function of m Z . The curves are obtained by fitting the µ − τ couplings for the different Z masses, to account for ∆a µ and respecting the constraints from the ratios R Z αβ and R τ µe . We then fit the other couplings, g eτ L,R and g eµ L,R , in order to drive the value of the most stringent cLFV processes (i.e. µ → eγ and τ → µēµ) to their current experimental upper bounds (see Table 1).
These illustrative theoretical predictions (which do not take into account effects of operator mixing due to RG evolution for large m Z > Λ EW ) should be compared with the future sensitivity of the associated dedicated searches, see Table 1.
As can be seen on the left panel 11 of Fig. 10, and while the maximal predictions for the cLFV   tau decays 12 in general lie beyond experimental sensitivity (at best, BR(τ → eee)∼ O(10 −13 ) ), one has very good prospects for the observation of certain rare muon decays. This is the case of both µ → eγ, and muon-electron conversion in nuclei: the former maximal branching ratio is fixed to its current experimental upper limit (≤ 4.2 × 10 −13 ), and the maximal conversion rate for Aluminium nuclei ∼ O(10 −17,−16 ), potentially within future experimental sensitivity. The current experimental limit for BR(µ → eγ) also precludes sizeable rates for µ → eee (within future sensitivity), such that a future observation of both µ → eγ and µ → eee at comparable rates would falsify Z models under discussion 13 . Also notice that the (full) lines for µ → eγ and τ → µēµ appear superimposed with the corresponding dashed lines denoting the associated experimental bound; this is a consequence of maximising the couplings to be within 1σ of the current experimental bounds. For illustrative purposes, we also display the prospects for the maximal contributions to the cLFV Z decays; as can be seen, and for m Z ≈ 500 GeV, one has BR(Z → e ± τ ∓ )∼ 10 −10 , marginally at future FCC-ee sensitivity [137]. Although not included here, let us notice that the maximal expected rates for the other cLFV Z decays are BR(Z → µ ± τ ∓ )∼ 10 −19 and BR(Z → e ± µ ∓ )∼ 10 −16 .
On the right-handed panel of Fig. 10 we display the prospects of the maximal expected deviations in what concerns LFUV ratios of Z → decay widths, together with the SM predictions and the experimental data. It is important to notice that, with the exception of R Z τ µ , the Z-decay universality ratios are generically predicted to be slightly smaller than in the SM. This is due to sizeable interference effects between the Z -loop contribution and the SM tree-level diagrams.

A (light) flavour violating Z at a future muon collider
In addition to the vast array of indirect (and direct, low-energy) searches for a new Z boson, its presence has also been the object of dedicated programmes at high-energy colliders, from LEP [164] effect. 12 The maximal predictions for τ → µγ decays appear to be constant (i.e., independent of m Z ) since the relevant couplings are those also responsible for τ → µēµ. Since the latter rate is fixed at its current limit for each value of m Z , one recovers the observed behaviour. 13 Notice that the cLFV observables discussed, and in particular their falsifying role, generically probe the presence of a new leptophilic Z with flavour violating couplings (combinations of g µτ with either g µe or g eτ ), and not directly the potential of such constructions to offer an explanation to ∆aµ.
to current efforts at the LHC [165][166][167][168][169]. Likewise, extensive work had been done in what concerns the prospects for the discovery of such a NP mediator in the (near) future, in particular at (HL)-LHC [108,170], Belle II [162] and FCC-ee [108]. In recent years, an increasing interest for future muon colliders has been put forward by the high-energy particle physics community [111,114,171,172], and the prospects for the discovery of a Z in such facilities have been also considered (see [91]).
In general, lepton flavour violating Z bosons can be produced at colliders in processes such as ff → µ ± µ ± τ ∓ τ ∓ , with the striking final state signature of two same-sign lepton pairs [108,162,170]. In order to complement the discussion, we now address the possible smoking guns of a leptophilic flavour violating Z at muon colliders. In this section we thus explore possible off-resonance signatures arising from the presence of a cLFV Z at µ + µ − colliders, discussing both the associated production cross sections, as well as forward-backward asymmetries.
Here, we neglect the initial muon masses and therefore Higgs-exchange diagrams. Furthermore, effects of initial and final state radiation corrections are not taken into account as they are beyond the scope of the present study. For a detailed discussion of the importance of these effects see e.g. [173] and references therein.

cLFV Z production at a muon collider
In what follows, we consider distinct scenarios, both concerning the operating centre of mass energy (i.e. √ s) and m Z . In each case, the couplings are determined as to comply with the best-fit regimes for the flavoured observables so far considered: saturating ∆a µ , and complying with all the relevant cLFV and LFUV bounds. In view of the discussion of the previous sections, the very stringent associated constraints lead to tiny Z − e − couplings (with = µ, τ ); we thus focus on tau-pair production at a future muon collider. In the present model, the process µ + µ − → τ + τ − receives contributions from SM γ and Z boson exchanges (s-channel) as well as from t-channel Z exchange, leading to the following matrix element: in which we highlight the sign difference between the SM (s-channel) contributions and the Z t-channel one. The differential cross section is given by with M defined in Eq. (58), and the details of the computation being collected in Appendix A.3. In the above, θ is the final state lepton angle with respect to the colliding muon direction in the di-tau centre of mass frame. Concerning the phase space integration, we employ a (realistic) angular cut of −0.99 ≤ cos θ ≤ 0.99 to account for the finite detector volume. We stress here that θ does not correspond to a reconstructed angle from a detector simulation and τ identification. The results concerning the production cross section σ(µ + µ − → τ + τ − ) as a function of √ s are presented in Fig. 11, for different values of m Z . The SM predictions are separately displayed -we notice that the computation of the latter is done only at leading order and that we neglect the muon mass. In all cases, the values of the couplings g µτ L,R are determined by fitting them to the likelihood containing (g − 2) µ , the Z-decay LFU ratios R Z αβ and the τ -decay LFU ratio R τ µe . The results of these fits, for selected values of m Z are shown in Table 2. Subsequently, for a fixed mass m Z and √ s, we sample the couplings around the best fit point according to the likelihood, in order to estimate the uncertainty of our predictions for µ + µ − → τ + τ − processes. In the following, the shaded regions denote the uncertainty corridor corresponding to the interval between the 16 th and 84 th quantiles of the samples. While for light Z mediators (such as those considered in most of the numerical analysis of the previous section) the behaviour of the current BSM construction hardly deviates from the SM expectation, the situation becomes visibly different for heavier states. This is a direct consequence of fitting the g µτ L,R couplings (which thus vary for every value of m Z ), and these -as seen from Fig. 9 - significantly increase for m Z 30 GeV, especially the right-handed couplings g µτ R (see Table 2 below). The latter also present the largest uncertainty in the fit, especially for m Z 150 GeV, as visible from Fig. 9, and this accounts for the spread in the different (shaded) regions associated with the coloured curves.  Table 2: Best fit values and associated uncertainties of the µ − τ couplings for different Z masses; the fit includes (g − 2) µ , R Z αβ and R τ µe constraints.
Studies of heavy Z resonances (m Z ∼ O(TeV)) at the LHC have been recently pursued, also emphasising the constraining role of the forward-backward asymmetry [174]; however, due to s-channel exchange, the studied angular distributions exhibit a resonant behaviour at the pole of the new boson mass, while deviations for lower √ s, e.g. at around the Z-pole, are negligible. In contrast, and as we will subsequently discuss in the present study, deviations in the forward-backward asymmetry A FB are quite significant already at lower energies, and considerably depart from the SM predictions at large √ s, due to the t-channel exchange of the new boson. The forward-backward (muon) asymmetry has also been explored in context of lepton flavour violating (eµ) axion-like particles aiming at explaining ∆a µ ; these light mediators are also responsible for new t-channel contributions in e + e − → µ + µ − scattering, searched for at the Belle II experiment [175].
Complementary insight can be obtained by considering the dependency of the σ(µ + µ − → τ + τ − ) production cross section on m Z , for different values of the centre of mass energy, √ s. This is displayed on the panels of Fig. 12, in which we again compare the SM predictions with those of the Z extension under consideration. As manifest from Fig. 12, for sufficiently "heavy" Z bosons (i.e. m Z 30 GeV) one has a clear departure from the SM prediction; for most choices of √ s, and as already seen in Fig. 11, , rendering this observable highly sensitive to the presence of heavier Z bosons; even with low statistics, a strong signal in conflict with the SM prediction could be expected. This information is summarised in Table 3, in which we collect the predictions for the µ + µ − → τ + τ − production cross section 14 for several values of m Z and √ s.  The benchmark values for √ s = 0.5, 1, 3, 10 TeV have been chosen in agreement with recent muon collider design proposals [112,114,171,172]. With the anticipated integrated luminosity of O(ab −1 ) [114,171,172], the µ + µ − → τ + τ − cross section can be expected to be precisely measured at such a facility.
A more extensive study making use of the usual Monte-Carlo "tool-chain" and detector simulation is left for subsequent studies.

Forward-backward asymmetry A FB
The forward-backward asymmetry, A FB , has been extensively used a powerful probe for the presence of NP in association with numerous production modes. At LEP-1 and LEP-2, scans of √ s around the Z-pole have been performed, in which different final states were studied [176]. Most of the results have been found to be in excellent agreement with the SM predictions. Relying on several simplifying assumptions, we have evaluated how the cLFV Z boson under study could lead to a departure from the SM predictions concerning the forward-backward asymmetry associated with the µ + µ − → τ + τ − process at a future muon collider. We re-iterate that ours is strictly a simple (naïve) tree-level calculation, without taking into account potentially relevant initial and final state radiation (ISR and FSR). Likewise, the results for A FB around the Z-pole mass are displayed strictly for illustrative purposes. The forward-backward asymmetry is defined as with in which the upper (lower) integral limits, "max (min)", denote the extreme values of cos θ considered. The prospects for A FB (µ + µ − → τ + τ − ) as a function of the centre of mass energy are presented in Fig. 13, for several values of m Z , and compared to the SM expectation. As can be seen, there is a clear difference between the SM and the Z predictions, increasing with larger √ s. Leading to these plots, we have followed the same procedure as in the previous subsection -considering for every value of m Z the best-fit intervals for the couplings g µτ L,R (to comply with cLFV and LFUV bounds, and saturate ∆a µ ). As previously mentioned, in our study we have also imposed (realistic) cuts on the integration limits (which also indirectly reflect finite detector effects), see Eq. (61): we have taken −0.99 ≤ cos θ ≤ 0.99. Although A FB (µ + µ − → τ + τ − ) is significantly larger (and distinguishable from the SM expectation) for √ s 200 GeV, the predictions for different m Z values become indistinguishable for √ s ∼ 1 TeV; also recall that, as shown in Fig. 11, the Z -induced flavour violating rates become very small for high energies (as a consequence of the associated t-channel event topology). Further -more detailed -computations might be required in this case.
It is also interesting to notice that there is a clear deviation from the SM expectation for √ s below the Z-boson pole: this is manifest from inspection of the right panel of Fig. 13, where one can easily have sizeable deviations in such a √ s regime.
Similarly to what was done in the previous subsection, in Fig. 14  Notice that for intermediate regimes of √ s, the expected A FB (µ + µ − → τ + τ − ) comes closer to its SM prediction for large m Z . This can be understood from the ∆a µ driven fit, due to which the couplings do not grow proportionally to the mediator's mass. Thus, for large m Z and comparatively "small" √ s, the new contributions are sufficiently suppressed. The striking differences with respect to the SM expectation are summarised in Fig. 15, in which we compare the prospects for the production cross section and A FB , considering As can be seen, our results show that even for comparatively small √ s the relevant regions of the model's parameter space can be falsified as an explanation to ∆a µ at a future muon collider. At larger √ s, precise measurements of the µ + µ − → τ + τ − production cross section may be used to constrain the mediator mass, even in the absence of a resonance in the spectrum.

Further discussion and outlook
Motivated by a minimal framework to address the tensions between theory and experiment in the anomalous magnetic moments of (light) charged leptons, we have considered a SM extension via a Z model. These well-motivated constructions have been extensively investigated in the literature (also in view of their potential interest in relation with the so-called LFUV anomalies in B-meson decays), and their flavour conserving couplings to fermions are very strongly constrained by a vast array of experimental measurements and searches. In our study, we have thus considered an ad-hoc (toy) model of a leptophilic light Z , which only couples in a flavour-violating manner to leptons. As we have argued here, ∆a Rb e and ∆a Cs e can be separately accounted for (for example, respectively relying on a same-sign or opposite-sign combination of the g eτ L,R couplings); however, a simultaneous explanation of ∆a e (either ∆a Rb e or ∆a Cs e ) and ∆a µ is precluded due to having excessive contributions to numerous cLFV (and also LFUV) observables. In particular, the most stringent bounds emerge from cLFV τ → µēµ decays (for g eµ L,R ), and from µ → eγ decays (g eτ L,R ); for g µτ L,R couplings, the LFU-probing ratios, R Z αβ (α = β) and R τ eµ play the most constraining roles.  Thus, and focusing on an explanation to the current tension in the muon anomalous magnetic moment, one is led to ∆a µ -favoured regimes, both for the relevant couplings and for the mass of the new mediator. In all cases, the latter regimes -and the strong constraints from cLFV and LFUV observables -strongly suggest that in this class of models one expects SM-compatible values for ∆a e . For the tau magnetic moment, one predicts values compatible with 0, typically ∆a τ ∼ [−10 −6 , −10 −7 ].
In what concerns cLFV observables, and in spite of their very constraining role, one has good prospects for the observation of certain rare muon decays in the regimes favoured by an explanation of ∆a µ , in particular for µ → eγ, and muon-electron conversion in nuclei. Interestingly, in these ∆a µ -favoured regimes, the stringent constraint imposed by the current experimental limit on BR(µ → eγ) leads to a theoretical upper limit for BR(µ → eee) that is smaller than the sensitivities of next-generation experiments. Thus, a future observation of the latter two cLFV observables would contribute to falsify these Z models.
We have also emphasised the role of a future muon collider in probing the model's regimes which are preferred by ∆a µ , for a wide range of Z masses. Relying on a simple estimation of the production cross-section, we have pointed out that one should have sizeable σ(µ + µ − → τ + τ − ), for various choices of the centre of mass energy. In addition to deviations from the SM regarding σ(µ + µ − → τ + τ − ), the forward-backward asymmetry A FB (µ + µ − → τ + τ − ) is also expected to exhibit distinctive features, which can be used to put this minimal construction to the test.
In the future, a more precise determination of a exp µ , possible in view of the expected reduction in the FNAL measurement uncertainty (ten-fold improvement in statistics [177]), together with improved SM predictions, will further allow to constrain the present model's parameter space. Likewise, clarifying the situation in what concerns the electron (g − 2) would also contribute to further test (or falsify) the model.
The toy model here considered -and despite not succeeding in providing a simultaneous explanation to the apparent tensions in ∆a e,µ -offers an interesting starting point to the analysis of complete SM extensions via new U (1) mediators. As pointed out, new fields should be added -as is the case of a new scalar responsible for symmetry breaking and for giving a mass to Z boson, and additional (heavy) fermions to properly cancel gauge anomalies and provide a realistic origin of the (strictly) flavour violating couplings. In turn, one expects that new interactions of the scalar with SM matter (i.e. the new Yukawa couplings) will then contribute to the many processes described above; one must then explore if in this case it is possible to account for both (g − 2) e,µ without conflict with the relevant cLFV and LFUV observables which so far preclude an simultaneously explanation of the tensions. Further extensions are also possible concerning the fermion (matter) sector, such as introducing right-handed neutrinos, and considering potential neutrino mass generation mechanisms. Likewise, one can also envisage considering the role of the new fields and couplings in what concerns an explanation to the dark matter problem.

Acknowledgements
This project has received support from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No. 860881 (HIDDeν network) and from the IN2P3 (CNRS) Master Project, "Flavour probes: lepton sector and beyond" (16-PH-169).

A Appendix
A.1 New contributions to LFC Z → + − decays In this appendix we summarise the most relevant details concerning the new contributions to the Z → + α − α decay (same-flavour charged lepton pair), arising from a Z in the loop. The matrix element can be written as where u α , v α denote the α spinors, µ (q) the Z-boson polarisation vector with q the Z momentum, and C V,T X are the vector and tensor coefficients (computed from the decay triangle diagram). For Hermitian (real in our case) Z couplings, these are given by: In the above, denote the Passarino-Veltman functions (with P V = 0, 1, 2, 00, 11,12,22) and m i the mass of the internal lepton. Notice that the dimensions of the vector and tensor coefficients are respectively C V L = m 0 and C T L = m −1 (recall that certain Passarino-Veltman functions are dimensionful quantities). Finally, the RH coefficients can be obtained by exchanging L ↔ R in the expressions for C V,T L . We further set q 2 = m 2 Z for the on-shell decays. Moreover, and as already discussed in the main part of this work, we treat the divergences appearing in the B 0 and C 00 functions in the MS-scheme and set the 't Hooft scale to µ 2 = m 2 Z . The contributions to the leptonic Z decay width (SM-like, Z and interference terms) can be cast as with and Γ SM full is given in [125] at 2-loop accuracy. In the above, the Källén-function is defined as

A.2 Lepton flavour violating Z decay amplitude and coefficients
We now consider the effects arising from the presence of the new Z boson in what concerns cLFV Z → + α − β decays. Analogously to the previous appendix, the matrix element can be cast as: already explicitly cast in terms of the vector and tensor coefficients, C V,T L . In the above, u β , v α denote the β , α spinors respectively, µ (q) is the Z-boson polarisation vector, with q the Z momentum. As before, K V,T X can be computed from the decay triangle diagram and are given, for hermitian (real) Z couplings, by: ) +g iα, * L g iβ L g Z L (−(m 2 β (m 2 i − 2m 2 Z )) + m 2 α (2m 2 β − m 2 i + 6m 2 Z )) C 1 + 4g iα, * R g iβ R g Z R m α m β m 2 Z + g iβ R g iα, * L m β m i (−2g Z R m 2 Z + g Z L (−m 2 α + m 2 i − 2m 2 Z )) +g iα, * R g iβ L m α m i (−2g Z R m 2 Z + g Z L (−m 2 β + m 2 i − 2m 2 Z )) +g iα, * L g iβ L g Z L (−(m 2 β (m 2 i − 6m 2 Z )) + m 2 α (2m 2 β − m 2 i + 2m 2 Z )) C 2 + m α −(g iβ R g iα, * L (g Z R + g Z L )m α m β m i ) − g iα, * R g iβ L (g Z R m 2 α + g Z L m 2 β )m i +g iα, * R g iβ R m β (g Z L m 2 i + g Z R (m 2 α + 2m 2 Z )) + g iα, * L g iβ L m α (g Z R m 2 i + g Z L (m 2 β + 2m 2 Z )) C 11 As in the previous appendix, C P V ≡ C P V (m 2 α , q 2 , m 2 β , m 2 Z , m 2 i , m 2 i ), are the Passarino-Veltman functions, i is the internal lepton, and the RH coefficients are obtained by exchanging (L ↔ R). As before, the dimensions of the vector and tensor coefficients are respectively K V L = m 0 and K T L = m −1 . As before, we treat the divergences appearing in the B 0,1 and C 00 functions in the MS-scheme and set the 't Hooft scale to µ 2 = m 2 Z .
Finally, the partial width for the decay Z → + α − β reads: in which "c.c." denotes the complex conjugate.
A.3 Tau-pair production at a muon collider: matrix elements for µ − µ + → τ − τ + The matrix element for the µ − µ + → τ − τ + process, considering s-channel Z and γ exchange, as well as t-channel Z contributions (and neglecting the muon mass) is given by as introduced in Eq. (58), with The matrix elements squared (and interference terms) are accordingly given by: where e is the electric charge, g cw g A(V ) are the axial (vector) couplings of the Z boson to charged leptons, Γ V denotes the width of the vector boson, with V = Z, Z (see Appendix A.4 for details). The Mandelstam variables, s, t, u, can be expressed in function of the energy √ s and the angle θ (the final state lepton angle with respect to the colliding muon direction in the di-tau centre of mass frame) as:

A.4 Width of the Z
The Z boson width can be estimated from its tree-level decays to charged leptons and neutrinos. The decay rate of a vector boson V to two fermions F 1,2 is given by (see, for example [178]): (85) in which λ(a, b, c) is the Källén-function already defined in Eq. (69). For Z → α β decays, the couplings b L,R are given by b L,R = g αβ L,R , whereas for Z → ν α ν β decays one has b L = g αβ L (with b R = 0). The total width Γ Z is then determined by the sum over all possible leptonic final states, neutral and charged, flavour conserving and/or flavour violating.