Gauging the accidental symmetries of the Standard Model, and implications for the flavour anomalies

We explore the possibility that lepton family numbers and baryon number are such good symmetries of Nature because they are the global remnant of a spontaneously broken gauge symmetry. An almost arbitrary linear combination of these symmetries (together with a component of global hypercharge) can be consistently gauged, if the Standard Model (SM) fermion content is augmented by three chiral SM singlet states. Within this framework of $U(1)$ extensions of the SM one generically expects flavour non-universality to emerge in the charged leptons, in such a way that naturally prevents lepton flavour violation, by aligning the mass and weak eigenbases. For quarks, all the SM Yukawa couplings responsible for their observed masses and mixings arise at the renormalisable level. We perform fits to show that models in this class can explain $R_{K^{(\ast)}}$ and the other neutral current $B$ anomaly data if we introduce a heavy vector-like quark to mediate the required quark flavour violation, while simultaneously satisfying other constraints from direct $Z^\prime$ searches at the LHC, $B_s$ meson mixing, a number of electroweak precision observables, and neutrino trident production. Within this symmetry-motivated framework of models, we find interesting implications for the flavour anomalies; notably, any axial couplings of the $Z^\prime$ to electrons and muons must be flavour universal, with the flavour universality violation arising solely from the vector-like couplings. We also comment on the generation of neutrino masses in these models.


Introduction
The renormalisable Standard Model (SM) lagrangian possesses a number of accidental continuous global symmetries, namely baryon number symmetry, U (1) B , and three individual lepton number symmetries, U (1) e , U (1) µ , and U (1) τ . These accidental symmetries are a great success of the SM, since they appear to be extremely good symmetries of Nature that are borne out by (almost) all particle physics experiments to date. The bounds on baryon number-violating processes, most famously from proton decay [1], and on lepton flavour violating (LFV) processes such as µ → eγ, µ → eee, and τ → µµµ, are very strong [1], and are forecast to strengthen further at future experiments [2]. There is only one hint from the world of particle physics to suggest that any of these accidental symmetries are in fact broken, and that is the inference of non-zero neutrino masses from the observation of neutrino oscillations. 1 This therefore begs the question: why are these such good symmetries of Nature? As soon as we begin to write down higher-dimension operators in the SM effective field theory (SMEFT), we find that these accidental symmetries of the renormalisable SM lagrangian are broken. The lepton number symmetries are broken already at dimension five, by the inclusion of a Weinberg operator of the form (LH) 2 /Λ, where Λ denotes the cut-off scale of the SMEFT. Baryon number is broken at dimension six by four-fermion operators such as QQQL/Λ 2 .
Thus, the first possibility is that U (1) B and U (1) L i are such good symmetries simply because the cut-off Λ of the SMEFT is very high with respect to energy scales E probed by current particle physics experiments, such as the LHC, so that even dimension five operators are suppressed by a factor E/Λ 1. While such a scenario certainly remains a viable possibility, it is a rather pessimistic one for particle physicists to swallow, for it implies, by and large, that any new physics beyond the Standard Model (BSM) is likely lying out of reach of current colliders.
The second and more intriguing possibility is that, in some extended BSM theory that arises due to new physics at a scale Λ LFUV , these global symmetries remain accidental at the level of the renormalisable lagrangian. If this is the case, one would expect a natural separation between the scale Λ LFUV of those higher-dimension operators which respect these global symmetries, and the scale Λ LFV Λ LFUV of those operators which violate them (where the reason for these names shall soon be apparent). 2 In such a scenario, the true cut-off Λ LFUV of the SMEFT can in fact be brought much lower in energy than the scale at which U (1) B and U (1) L i are violated. While this second option is certainly favoured by the biases of optimistic model builders, until recently there has been no physics case for favouring either one of these two hypotheses over the other.
Fortunately, this impasse is now being challenged by the recent emergence of a set of intriguing discrepancies between SM predictions and experimental measurements, in the neutral current decays of B mesons, which would seem to favour the second of the two possibilities just discussed. For example, the ratio of branching ratios R K ( * ) ≡ BR(B → K ( * ) µ + µ − )/BR(B → K ( * ) e + e − ) is equal to unity in the SM to the percent level, for di-lepton invariant mass squared bin q 2 ∈ [1.1, 6] GeV 2 , but LHCb has measured [3,4]  understand the matter-antimatter asymmetry, and thus the evolution of structure in our Universe. 2 One might try to pursue a middle way between these two options, and suggest that there exists such a hierarchy of scales in the SMEFT, with U (1)B-or U (1)L i -violating operators being suppressed by some cut-off scale ΛLFV ΛLFUV, not because of any underlying gauge symmetry, but simply because of small couplings in the EFT (in other words, one chooses to forgo the expectation from naïve dimensional analysis (NDA), and the 'second scale' really signifies the existence of very small couplings c such that the ratio ΛLFUV/c ∼ ΛLFV ΛLFUV). This approach is at least consistent, in the sense that the smallness of the couplings c can be radiatively stable in the low-energy EFT. R K * = 0.69 +0. 11 −0.07 ± 0.05 in this q 2 bin, corresponding to deviations from the SM by approximately 2.5σ each. LHCb has also measured R K * = 0.66 +0. 11 −0.07 ± 0.03 for the low momentum bin q 2 ∈ [0.045, 1.1] GeV 2 , which is again about 2.5σ under the SM prediction [3]. There are further notable discrepancies with the SM predictions in measurements of BR(B s → µµ) [5][6][7][8], and in B → K * µ + µ − angular observables such as P 5 [9][10][11][12]. For a comprehensive survey of these anomalies in the decays of neutral B mesons, which we henceforth refer to collectively as the 'neutral current B anomalies' (NCBAs), see e.g. [13].
Taken together, these measurements point towards lepton flavour universality violation (LFUV) between e and µ. 3 But crucially, there is no evidence for LFV. Thus, these tantalising first hints of BSM physics at the LHC seem to respect the accidental global symmetries of the SM, despite violating lepton flavour universality. The magnitude of the deviations from the SM in the various NCBAs suggests there is new physics at a cut-off scale of order Λ LFUV ∼ 30 TeV [13,[21][22][23][24][25] or thereabouts (most pessimistically, constraints from perturbative unitarity imply the new physics scale cannot be larger than about Λ LFUV ∼ 80 TeV [26]). This energy scale is probed at the LHC by measurements of many rare flavour observables, not just those mentioned above, and there are no other signs of new physics at this scale, and certainly no evidence for the violation of U (1) B or U (1) L i . In this way, if the NCBAs persist in the wake of future measurements (from LHCb, Belle II, and others [27]), and if these or any other new anomalies continue to respect U (1) B and U (1) L i , it would seem to suggest that the SM's global symmetries U (1) B and U (1) L i are likely a consequence of some underlying dynamics, such as a gauge symmetry.
If we adopt such an optimistic interpretation of the NCBAs, then which symmetry might one gauge, which can both explain the NCBAs and would also underwrite these global symmetries? At the level of the SMEFT, the NCBA data can be explained by BSM contributions to the following four-fermion operators where the SM contribution is C e,SM L = C µ,SM L 8.64/(36 TeV) 2 and C e,SM R C µ,SM R = −0.18/(36 TeV) 2 (borrowing the numerics from Ref. [28]), which is due to one-loop W exchange. Global fits to the data favour a significant BSM contribution to C µ L , and are consistent with non-zero BSM contributions to the other three operators [13,[21][22][23][24][25]. One simple possibility is that all these operators receive BSM contributions due to the tree-level exchange of a heavy Z vector boson, which couples to muons and (possibly) electrons, in addition to possessing a flavour violating coupling to bs. While the NCBAs can also be explained by tree-level leptoquark exchange or by various loop-induced processes, in this paper we shall restrict our attention to Z models, in which the Z arises from a spontaneously broken, flavour-dependent U (1) gauge symmetry, which we shall denote henceforth by U (1) X .
There are many possible U (1) X symmetry groups that can be gauged in order to explain the NCBAs, and the literature on such models has grown vast (see e.g. Refs.  for an incomplete list). In this work, we shall make a number of simplifying assumptions before arriving at a new framework of SM×U (1) X theories in which the protection of the SM's accidental symmetries is built-in. Firstly, we explore only U (1) X charge assignments which allow a SM-like Yukawa sector for the quarks, at the renormalisable level. Within such an approach we make no attempt to shed light upon the peculiar flavour structure of the SM, which features large hierarchies in fermion mass parameters and the quark mixing angles. On the other hand, we avoid having to delve more deeply into the model building to explain the origin of Yukawa couplings if they are to be banned at the renormalisable level, as they are in e.g. Refs. [53][54][55].
In the lepton sector, we shall require that only the diagonal charged lepton Yukawa couplings are U (1) X -invariant, and thus present in the renormalisable lagrangian. This is achieved (given only a single Higgs, possibly charged under U (1) X ) by requiring the different family leptons have different U (1) X charges, but that the differences between the left-handed and right-handed lepton charges are family universal. Such a charge assignment automatically implies LFUV, since the Z couples differently to each lepton family. Furthermore, in banning all the off-diagonal charged lepton Yukawa couplings (at the renormalisable level), the charged lepton mass eigenbasis is aligned with the weak eigenbasis, 4 and hence no lepton flavour violating neutral currents will be induced in the physical mass basis. Thus, the very same U (1) X gauge symmetry which introduces LFUV also prevents LFV, in opposition to the claims made in Ref. [56]. 5 We must also require our U (1) X symmetry be free of all gauge anomalies, including mixed and gravitational anomalies. 6 In these senses, the SM×U (1) X theories we shall define may be regarded as technically 'complete', minimal extensions of the SM, albeit with an unexplained flavour structure as in the SM itself.
What is the most general SM×U (1) X theory which satisfies these criteria? In the spirit of minimality, we want to add as few chiral states as possible beyond those of the SM. The most minimal option, of course, is to add no BSM states at all (beyond the Z , and a heavy scalar field Φ introduced to spontaneously break U (1) X ). In this case, it has been known for a long time that the only anomaly-free U (1) X charge assignments consistent with our criteria correspond to gauging L i − L j , the difference of any pair of lepton numbers [59]. 7 there is not a global part of that gauge symmetry that remains unbroken (when considering only its action on fermions). This is the case, for example, with the global part of the hypercharge symmetry of the SM after electroweak symmetry breaking. 6 We note in passing that in many Z models which are necessarily only low-energy EFTs (for example, because they do not permit a gauge-invariant Yukawa sector), it is not essential to cancel gauge anomalies at low energies; it might be possible in such models for anomaly cancellation to be restored in the high-energy theory by e.g. 'integrating in' a set of heavy chiral fermions, or by the Green Schwarz mechanism [57,58]. In our setup, we seek to embed our Z in a fully renormalisable extension of the SM, and thus anomaly cancellation is for us an essential requirement. 7 In fact, while this statement is 'common lore' amongst physicsts, we shall discover in §2 that, under these The particular choice of gauging L µ − L τ offers a compelling Z explanation of the NCBAs (after flavour violating couplings to quarks are introduced, for example through effective non-renormalisable interactions), which moreover underwrites the SM's global lepton number symmetries with a gauge symmetry in the manner described above [31,37]. The Z boson of such a model can even mediate interactions with a dark sector; connecting the NCBAs with the dark matter problem through such a gauged L µ − L τ leads to a highly-testable model with rich phenomenology [60]. The next minimal option is to introduce a 'dark sector' of SM singlet chiral fermion states, which are charged only under U (1) X . We shall here consider augmenting the SM matter content with up to three such states. Since the extra states are chiral only with respect to the U (1) X symmetry, they do not spoil the cancellation of anomalies in the SM gauge sector. Nonetheless, these states contribute to two out of the six anomaly coefficients involving U (1) X , and their inclusion can thus be used to cancel two anomaly coefficients that would otherwise be non-zero, thereby opening up a wider space of anomaly-free charge assignments beyond just L i − L j . Within this setup, there is a four-parameter family of such anomaly-free U (1) X symmetries that can be gauged, generated by where a e , a µ , a τ and a y are rational coefficients, and T Y denotes the generator of (global) hypercharge. We may rewrite this generator in terms of the generators of the accidental global symmetries of the SM, as where here a i ∈ {a e , a µ , a τ }. This intriguing result, which we shall review in §2, was originally derived in Ref. [61], in which the three chiral dark states were interpreted as right-handed neutrinos. 8 To summarise, if one allows the addition of three SM singlet states to 'soak up' anomalies, then the most general anomaly-free U (1) X charge assignment for the SM fermions, which allows a fully generic quark Yukawa sector and a strictly diagonal charged lepton Yukawa matrix, corresponds to gauging an almost arbitrary linear combination of the (otherwise accidental) global symmetries of the SM (including the 'global part' of hypercharge symmetry). The word 'almost' indicates an important caveat, that the linear combination in (1.3) is not entirely arbitrary; in particular, the component of baryon number is fixed by conditions, one may gauge a slightly larger symmetry, of the form (Li −Lj)+aY , where Y denotes hypercharge and a is any rational number. 8 A three-parameter subset of these gauge symmetries, in which it was assumed that aY = 0, was also considered in Ref. [62] (and subsequently in Refs. [63,64], for example), with the goal of explaining the neutrino masses and mixing parameters. In our setup, we shall not interpret the three chiral dark states as right-handed neutrinos; rather, we shall suggest that neutrino masses may arise from a rather different mechanism -see Appendix A. Various other models in which more than one of the global SM symmetries is gauged have also been considered in the literature, e.g. in Refs. [65,66].
the components of each of the lepton numbers. This means, for example, that one cannot gauge B + L, which is of course anomaly-full. As long as the coefficients {a e , a µ , a τ } are all different, gauging such a symmetry implies LFUV, while preventing LFV.
In the rest of this paper, we develop Z models based on gauging U (1) X symmetries in the family defined by (1.2), and explore their phenomenology. In some sense, these models provide a wide generalisation of the L µ − L τ model, and so several aspects of our setup shall be borrowed from Ref. [31] (such as the mechanism for generating quark flavour violating couplings of the Z to bs). We show in §4 that these models make a distinctive prediction for the structure of the NCBAs, which may be decomposed into a flavour non-universal vector coupling to leptons, plus a flavour-universal axial component (if a Y = 0). Furthermore, if this axial component of the anomaly is non-vanishing, the Z inevitably must acquire tree-level couplings to valence quarks. We extract new global fits to the NCBA data using flavio [67], in terms of the parameters of our model. We identify a physically motivated benchmark scenario in order to interpret these multi-dimensional fits, by fixing a µ = 1 and a τ = 0, and compute other important phenomenological bounds in §5 at this benchmark point in our parameter space, before concluding.
Our goal in this paper is not to explore the phenomenology of these theories completely; rather, we are content to demonstrate that our framework leads to rich phenomenology, as hinted at above, and that there is interesting allowed parameter space, thereby laying the groundwork for possible future studies.
In Appendix A we discuss how light neutrino masses might naturally arise within our framework, which induce lepton flavour violation at the higher scale Λ LFV . This involves a more in-depth examination of the dark sector of the theory, for the more enthusiastic of our readers.

A new framework for LFUV without LFV
We consider an extension of the SM by a flavour-dependent U (1) X gauge symmetry, under which all SM fields plus three SM singlet chiral fermions (which we will denote by ν i R ) may be charged a priori. The U (1) X symmetry will be spontaneously broken by the vacuum expectation value (vev) v Φ of a SM singlet scalar field Φ, at around the TeV scale, leading to a heavy Z gauge boson. We require: 1. That the assignment of U (1) X charges is anomaly-free.
2. That all Yukawa couplings in the quark sector are permitted at the renormalisable level.
3. That the flavour-dependent U (1) X gauge symmetry protects each individual lepton number symmetry, and thus forbids lepton flavour violating (LFV) processes. This is achieved by aligning the mass eigenbasis of charged leptons with the weak eigenbasis. Thus, we require that only the diagonal elements are present in the charged lepton Yukawa matrix (at the renormalisable level), which in turn requires lepton flavour universality violation (LFUV) between all three families.
We write the fermion fields (including the three SM singlets) as the following representations of the gauge group SU where the index i ∈ {1, 2, 3} denotes the family. Our assumptions regarding fully generic quark Yukawa couplings imply the quark charges under U (1) X are flavour universal, whereas our mechanism for preventing LFV requires the lepton charges are flavour non-universal. The scalar sector of the theory contains the Higgs and the scalar Φ, which carry the representations The charges under U (1) X , denoted by the set ofQs, are all assumed to be rational numbers. 9 In such an extension of the SM by a (flavour-dependent) U (1) X gauge symmetry, there are six independent anomaly coefficients which must vanish, corresponding to the six (potentially non-vanishing) triangle diagrams involving at least one U (1) X gauge boson: , and the mixed gauge-gravity anomaly involving U (1) X . Thus, anomaly cancellation implies a system of six non-linear equations over the eighteen chiral fermion charges listed above. 10 By rescaling the gauge coupling, we may take these rational charges to be integers, and then apply arithmetic methods (such as those introduced in Ref. [68]) to solve the resulting system of non-linear Diophantine equations.
In fact, it shall turn out that because of the heavy restrictions we are enforcing from the Yukawa sector, the subspace of anomaly-free solutions that we are interested in is essentially picked out by a strictly linear system of equations, and so may be extracted using basic linear algebra; thus, employing the Diophantine methods outlined in Ref. [68] in this case would be something like overkill.
We begin by enforcing the constraints on the charges coming from the Yukawa sector. As already noted, requiring a renormalisable quark Yukawa sector implies that q L , u R , and d R have flavour-universal U (1) X charges, which we denote byQ q ,Q u , andQ d . Then, requiring renormalisable Yukawas for up and down quarks, and diagonal charged lepton Yukawas in each generation, implies the following five linear constraints (which are all linearly independent) on the remaining ten charges:Q If the charges under U (1)X are not rational numbers, then such matter content cannot arise from any unified gauge theory with semi-simple gauge group G. There are good reasons to assume that such a unified theory ultimately describes the interactions of elementary particles in the ultraviolet. 10 One may of course add arbitrary scalars or vector-like fermions without affecting anomaly cancellation.
We are thus reduced, by these linear constraints, to a five-parameter family of charges. As might be expected, this family of charge assignments corresponds to gauging an arbitrary linear combination of the five accidental global symmetries of the SM Yukawa sector: U (1) B , U (1) e , U (1) µ , and U (1) τ , and (the global part of) hypercharge, U (1) Y . But such a charge assignment is not yet anomaly-free, an issue that we shall now remedy.
Firstly, it is helpful to notice that the three SM singlet dark states appear in just two of the six anomaly cancellation equations. Specifically, these are the (linear) gauge-gravity anomaly, whose coefficient is proportional to (2.2) and the (cubic) U (1) 3 X anomaly, with coefficient proportional to Thus, the charges of the SM fermions on their own must satisfy the other four anomaly equations independently.
Of these, three are linear, and moreover out of these three linear constraints, only one turns out to be linearly-independent from (2.1). 11 This additional linear constraint has the effect of fixing the component of baryon number in terms of the other symmetries (in such a way that excludes the gauging of B + L, which is well known to be anomaly-full). We are thus reduced by this subset of the anomaly cancellation equations to a four-parameter family of solutions. A particularly suggestive parametrisation for the SM fermion charges, in terms of four rational numbers a e , a µ , a τ , and a Y , is recorded in the first ten rows of the following: 11 This additional constraint can be taken to be either the U (1) 2 Y × U (1)X anomaly equation or the SU (2) 2 L × U (1)X anomaly equation or indeed a linear combination of these two, but not the SU (3) 2 × U (1)X anomaly equation.

Field
where in the bottom five rows we have also included the states in our 'dark sector' here for completeness; there are the three dark chiral fermion states ν i , which are SM singlets, the scalar Φ whose role it is to spontaneously break U (1) X by acquiring a vev, and a vectorlike heavy fermion denoted Q, which plays an important role in introducing quark flavour violation into the interactions of the Z , which we shall discuss in §3.1.
The coefficient of the U (1) 2 X × U (1) Y anomaly, which we have not yet considered, is proportional to the quadratic expression Somewhat surprisingly (or perhaps not), when we substitute into (2.5) the general solution parametrised in (2.4), we find that A quad = 0 for any values of (a e , a µ , a τ , a Y ) ∈ Q 4 .
At this stage, we turn to the two anomalies which are sensitive to the dark sector chiral fermions ν i R , via the anomaly coefficients in (2.2, 2.3). If we substitute in the charge assignment in (2.4), which we have derived as the most general charge assignment for the SM sector fermions (that is consistent with the other constraints from anomaly cancellation and the Yukawa sector), we obtain In the absence of the dark sector fermions (Q ν i =0), the only possible solution 12 to A cubic = 0 is to choose one of a e , a µ , and a τ to be zero, and the other two to sum to zero; a Y remains unconstrained. This also satisfies A grav = 0, and thus returns anomaly-free solutions of the form a(L i − L j ) + bY . But by introducing the three dark sector fermions, we can absorb any remaining anomalies in A grav or A cubic by ascribing them chargeŝ or any permutation thereof, which would correspond to the lepton numbers of right-handed neutrinos, with one assigned to each generation. 13 The anomaly-freedom of this charge assignment was originally shown in Ref. [61]. The charge assignment in (2.4) corresponds to gauging a U (1) X symmetry generated by i.e. an almost arbitrary linear combination of the accidental symmetries of the SM, namely baryon number, the three individual lepton numbers, and the global part of hypercharge (with that linear combination being 'orthogonal' to B +L). Of course, this is by no means a miracle; the very fact that these quantities are accidental symmetries of the SM lagrangian implies that the Yukawa sector is invariant under precisely these symmetries. What is surprising, at least to us, is that an (almost) arbitrary linear combination of these global symmetries can be made anomaly-free, and thus can be gauged, with only a minimal extension of the SM field content by three chiral SM singlets. Provided that a e , a µ , and a τ are all different, in other words that the charge assignment violates lepton flavour universality in all three families, the charged lepton Yukawa matrix will be strictly diagonal; thus, the same symmetry which introduces LFUV simultaneously prevents LFV.

Towards a model for the B anomalies
In this Section, we develop SM×U (1) X gauge theories in this family into phenomenological models capable of explaining the NCBAs. For ease of reference, we summarise the full 12 It is a rather famous theorem in number theory that a 3 e + a 3 µ + a 3 τ = 0, where ae, aµ, and aτ are three rational numbers, implies aeaµaτ = 0. 13 We are not claiming that the assignment of the dark sector fermion charges in (2.7) is the only solution to the pair of equations (2.6); rather, we are content that there is always guaranteed to be a solution to all the anomaly equations with dark sector charges of this form. Indeed, for certain rational values of the triple (ae, aµ, aτ ) it is known that other non-trivial solutions exist [69]. In particular, if we rescale the gauge coupling (without loss of generality) such that (ae, aµ, aτ ) and {Q ν i } are all integers, then an algorithmic method has been developed in Ref. [69] for finding all solutions to (2.6) for the six numbers (ae, aµ, aτ ,Q ν 1 ,Q ν 2 ,Q ν 3 ). Examples include (5, −4, −4, −1, −1, −1), (6, −5, −5, −3, −2, 1), and (11, −9, −9, −4, −4, 1) [69]. In the present paper we wish to be able to vary (ae, aµ, aτ ) freely within the rationals, in order to carry out fits to the NCBA data; in this situation, there do not necessarily exist any 'non-trivial' solutions to (2.6) beyond the charge assignment (2.7) that we choose. Nonetheless, we shall reconsider such 'non-trivial solutions' when we come to discuss neutrino masses in Appendix A. lagrangian of the model in Eq. (3.18). In order to mediate the flavour-changing neutral current interactions in (1.1), we must introduce flavour-changing quark couplings into our framework. Since the U (1) X gauge symmetry couples universally to the three quark generations, this cannot be achieved simply by the CKM mixing between the weak and mass eigenbases. The same issue afflicts (for example) models based on gauging L µ − L τ (in which the Z doesn't couple directly to quarks at all). We shall thus generate quark flavour violation in our model by following a similar procedure to that found in Ref. [31], which in turn is based on the 'effective operator' approach first suggested in Ref. [70]. We will then go on to describe the mass mixing between the neutral gauge bosons, specifically the Z and the Z , that occurs in these theories when a Y = 0, for which we largely follow Refs. [53,55].

Quark flavour violation
We introduce a heavy 14 vector-like quark field, denoted Q, whose left-and right-handed components both transform in the representation In other words, other than the U (1) X charge and its vector-like nature, the field Q is a 'heavy copy' of the quark doublet field q L in the SM. Indeed, it shall be convenient to decompose Q into its SU (2) L components, viz.
where the index anticipates that we shall soon view u 4 L and d 4 L as fourth family quark fields. Because Q is a vector-like fermion, it does not spoil the anomaly cancellation in our setup. Together with the fields already described above, this completes the field content of our framework of models, as recorded in (2.4).
Given the quantum numbers of Q and of the SM quark fields, we can write down the following terms in the lagrangian, which result in effective mass terms after both Φ and the Higgs acquire their vevs: The term −m Q QQ is simply a vector-like mass term for Q. The first two terms within the square brackets are the usual Yukawa couplings for the down and up quarks, while the second term in L mix leads to mass mixing between the vector-like quark and the SM quarks. Just like the SM Yukawa couplings (Y D ) ij , the Y Qi are (possibly complex) dimensionless Yukawa couplings, one for each down-type quark, which are parameters of the model. The U (1) X charge of Q R (and hence also Q L ) is then fixed by U (1) X -invariance of (3.2), to beQ Q = −Q Φ +Q q , as recorded in (2.4). To simplify our analysis, we shall assume the limit where m Q |v Φ Y Qi |. We may package together these terms into 4 by 4 quark mass matrices. We shall assume that we have first rotated in quark space to a basis where the SM 3 by 3 Yukawa matrix for the down quarks. The resulting 4 by 4 mass matrix for the down quarks is then given by where v is the usual Higgs vev), and the mixing parameters 1, given our limit of large m Q . We here use notation where the primed fields are current eigenstates, while the unprimed ones denote the physical mass eigenstates. We may diagonalise this 4 by 4 mass matrix completely with a further bi-unitary transformation of the form Terms of O( 2 ) in U L are in principle important when discussing flavour violation, but as we shall soon see, we do not in fact need to compute them. Concerning the rotation in the right-handed down quarks, we keep U R = 1, justified by the smallness of the quark masses compared to m Q . From these mixing matrices, we can compute, from now on neglecting higher order terms in the i parameters, the couplings of the Z 15 to down quarks within our model. Because the SM quarks have universal U (1) X charges, in the primed basis defined above (in which the SM Yukawa matrices Y D and Y U are already diagonalised), we have the following hadronic 15 In the following formulae, Z denotes the physical heavy gauge boson, which is a mass eigenstate. As we shall soon see, this is equal to the gauge field for U (1)X , which shall be denoted Xµ, up to small corrections couplings of the Z , where g X denotes the gauge coupling for U (1) X . Rotating to the mass basis, we obtain the following couplings of the Z to the three SM quark families, indexed as usual by i ∈ {1, 2, 3}, where the 3 by 3 matrices of Z couplings to the down quarks are given by The matrices of Z couplings to the up quarks are where V is the CKM matrix, and the index A ∈ {1, 2, 3, 4}.
To extract the flavour violating effects, which are of order O( 2 ), we can subtract from the 4 by 4 matrixL the flavour universal component equal toQ q δ AB . Using also the unitarity of U L , we then extract the flavour violating piece of L (d) ij to be (where i = j): where we used the fact thatQ Q −Q q = −Q Φ = −1. Of particular interest from the point of view of the NCBAs, there is a coupling of the Z to bs, of the form (3.11) Of course, within our setup there is no rotation between the mass and weak eigenbasis for the charged leptons because the U (1) X -invariant Yukawa terms are strictly diagonal. Thus, we only have diagonal couplings of the Z to charged leptons, with the charges as defined in (2.4), which result in LFUV without LFV.

Z − Z mixing
Provided a Y = 0, the Higgs carries U (1) X charge in our framework of models. This leads to mass mixing between the SM Z boson and the Z , which shall ultimately result in the Z inheriting some small flavour non-universality in its couplings to leptons. The universality of Z couplings to leptons is tightly constrained by precision measurements at LEP, which shall provide an important bound on our model when a Y = 0. To derive the following formulae for this Z − Z mixing, we closely follow Refs. [53,55].
As we have already set out, the U (1) X gauge symmetry is spontaneously broken by the SM singlet complex scalar Φ acquiring its vev. We shall here denote the original U (1) X gauge boson by X µ , reserving the name Z µ for the physical boson which is a mass eigenstate. The mass terms for the heavy gauge bosons come, of course, from the kinetic terms of the scalar fields H and Φ, where the covariant derivatives are Here, as usual, g and g denote the gauge couplings for SU (2) L and U (1) Y respectively.
Expanding the scalar fields about their vevs in (3.12), we find the following mass matrix for the neutral gauge bosons: (3.14) We rotate to the mass basis of physical neutral gauge bosons, where θ w is the Weinberg angle (such that tan θ w = g /g), and α z is the Z-Z mixing angle. The masses of the Z and Z boson are then the two non-zero eigenvalues of the above mass matrix. The mass of the Z is 16) and the mixing angle is where we are working in the limit that m Z m Z .
This concludes the description of our framework of models. In summary, the full lagrangian may be written as where L SM denotes the SM lagrangian (but with the Higgs kinetic and potential terms removed), X µν = ∂ µ X ν − ∂ ν X µ , 16 L mix may be read off from (3.2), L qZ from (3.6), L HΦ kin from (3.12), and L Z denotes the (flavour-diagonal, but flavour non-universal) Z coupling to leptons, given by where we have included here the Z couplings to the dark states ν i R . We summarise the Feynman rules associated with the important new physics couplings of the Z boson to SM fermions in Fig. 1, for the reader's convenience. Finally, V (H, Φ) is a scalar potential which determines the vevs of H and Φ (but which we do not specify beyond that), and L dark denotes any terms in the lagrangian involving dark states only. All that remains to specify the model fully is a discussion of how the neutrino mass sector may arise. Such a discussion involves a rather in-depth analysis of the dark sector of the theory (encoded in L dark ), and would be something of a digression at this stage, so we relegate this discussion to Appendix A.
The following two Sections are devoted to exploring the phenomenological consequences of this setup.

Implications for the B anomalies
After integrating out the heavy Z boson, we obtain BSM contributions to a host of dimension six Wilson operators in the SMEFT which are capable of explaining the observed NCBAs, depending on the values of the coefficients a e , a µ , a τ , and a Y . In §4.1 we will present the contributions to the relevant Wilson coefficients within our framework of models, which we shall see carry a particular structure allowing both vectorial and axial currents (but where the axial contributions must be lepton flavour universal). We then present results of our global fits to the NCBA data using flavio [67] in §4.2, after making a number of (physically well-motivated) simplifying assumptions to cut down our parameter space.

The anatomy of the B anomalies within our framework
At a generic point in the parameter space we are considering, the following four operators (all of which couple only to the left-handed bs current) are all present: 16 Note that we have assumed that any kinetic mixing between the Z and Z gauge fields is set to zero. where the coefficients are where α is the lepton flavour index α = {e, µ, τ }. We may convert these Wilson coefficients into the more conventional basis of vectorial and axial currents, corresponding to C α 9 ≡ (C α L + C α R )/2 and C α 10 ≡ −(C α L − C α R )/2, for which we find: This choice of basis makes clear that the NCBAs have an interesting, and simple, structure within our framework of models. A few noteworthy points are: • The LFUV must come entirely from the vectorial current.
• There may nonetheless be an axial current contribution, but this is lepton flavour universal. Indeed, this feature can be traced back to our initial assumptions regarding the renormalisable Yukawa sector, which implied that the differences between the lefthanded and right-handed lepton charges were family universal.
• The presence of this axial contribution requires a Y = 0, which then implies (i) that there is Z-Z mixing (sinceQ H = −a Y /2), and so important constraints from LEP lepton flavour universality measurements, and (ii) that there are necessarily couplings of the Z to valence quarks, as can be deduced from the charges in (2.4), and so important constraints from LHC direct searches in say pp → µµ.
• If we wish to remove couplings of the Z to quarks in order to loosen the constraints from direct searches, we require both that a Y = 0 (thus removing any axial component in the Wilson coefficients, as above) and that a e + a µ + a τ = 0. Thus, in this limit, there are only two independent parameters, which we can choose to be any two of (a e , a µ , a τ ).
• The expressions (4.4, 4.5) for the Wilson coefficients do not depend on the parameter v Φ , but only on the parameters from the gauge sector and from the mixing sector. The dependence on both v Φ and g X happens to cancel between the factors of M Z in the denominator and the couplings in the numerator.

Global fits at a benchmark point in the parameter space
Within this general class of models we have formulated, there is a large number of a priori free parameters. We have: • From the gauge sector: g X , a e , a µ , a τ , a Y .
• From the mixing sector: Y Qd , Y Qs , Y Qb , m Q .
• From the scalar sector: v Φ + . . . , where the . . . indicates additional parameters appearing the extended scalar potential, which we shall not be concerned with here. To thoroughly explore the phenomenology in all these parameters is a complicated task; in this paper we shall not attempt such a complete phenomenological characterisation of this parameter space. Rather, we prefer to make a number of well-motivated simplifying assumptions to cut down the parameter space. In this way, we shall define a "benchmark" region in our parameter space which is most relevant for the B anomalies and related phenomenology, and we shall find that there is room in this region to evade all current experimental constraints while remaining highly predictive.
The assumptions we make are as follows: 1. a µ = 1. This is in some sense a choice of normalisation, also made in Ref. [31], which forces a non-zero new physics contribution in the coupling of the Z to muons. While there remains a logical possibility for fitting the 'theoretically clean' ratios R K ( * ) with a µ = 0 and new physics only in the electron couplings, the inclusion of other NCBA data (for example P 5 in B → K * µµ decays, and BR(B s → µµ)) strongly favours new physics (or at least a sizable component) in the muon channel.
2. a τ = 0. The NCBAs concern only electrons and muons, and so are largely insensitive to the tauon couplings. Thus, the value of a τ shall have no effect on the global fits we will perform, and so it is convenient to set a τ = 0 at the outset.
3. Y Qd = 0. Fitting the NCBAs requires both Y Qs and Y Qb are non-zero, but does not require a non-zero coupling Y Qd . By choosing to set Y Qd = 0, we prevent new flavour violation beyond the SM in processes involving the down quark. With this choice we are therefore automatically safe from the bounds on e.g. kaon and B d mixing. Note that this assumption, like the others we are making, is not forced upon us by the data, but is a sensible choice which we make to reduce the parameter space.
4. g X is irrelevant for low energy processes. In the limit |Y Qi v Φ | < m Q which we are assuming, all the effects induced by the Z exchange are independent of g X (since two powers of g X from the gauge vertices in the numerator are always cancelled by two powers of g X in the denominator from the mass of the gauge boson m Z = g X v Φ ). If the mass of the Z is not too high, then the coupling g X and v Φ do indeed become two independent parameters, since bump searches for the Z depend on m Z itself. Nonetheless, if we are in the régime m Z m Z , 17 then direct searches for the Z only depend on the contact 4-fermion operators, independent of g X . We shall see in §5.1 that there is such a régime in our parameter space in which the Z couplings remain perturbative (in the sense that the Z width is less than, say, 30% of its mass, or thereabouts).

Y Qb Y *
Qs as a single parameter. Each of the couplings Y Qb and Y Qs appear only in the combination Y Qb Y * Qs in both the Wilson coefficients C α 9 and C α 10 given in (4.4, 4.5), and in the BSM contribution to B s mixing (see §5.2), which provides the main constraint sensitive to the quark flavour violation.
Given these assumptions, we perform global fits to the NCBA data in the plane of a e versus a Y , using flavio which, for the purpose of fitting, we write as The normalisation C is extracted point-by-point from the fit. The results of the fit are shown in Fig. 2. Details of the fit procedure are given in Appendix B. In the plot on the left-hand-side, only a subset of 'clean observables' (i.e. for which the theoretical uncertainties are smaller) are included in the fit. These clean observables are R K , R K * , and the branching ratios BR(B s → µµ), BR(B → X s µµ) and BR(B → X s ee). The dark (light) blue region extends to the 1σ (2σ) best fit contour. Interestingly, we find that a very large region of the parameter space in (a e , a Y ) allows a good fit to these clean LFU ratios. Indeed, we find two disconnected 'lobes' in the parameter space which fit the data at the 2σ level. In the plot on the right-hand-side, we further include all the branching ratios of the exclusive semileptonic branching ratios B → Kµµ, B → K * µµ, B s → φµµ and Λ b → Λµµ, as well as all CP -averaged angular observables in those decay modes into the fit. We see that the inclusion of additional variables beyond R K ( * ) has a dramatic effect on the fit, as indicated by the shaded 1σ and 2σ regions in green, where now a definite quasi-  Table 1. Best-fit values for a e and a Y , subject to the choices a µ = 1 and a τ = 0. The value C gives the overall normalisation of the Wilson coefficients at the best-fit point, see eqs. (4.6) and (4.7). The first line is for the fit to a subset of 'clean observables', corresponding to the left-hand-plot of Fig. 2.
The second line, for all observables, corresponds to the right-hand-plot of Fig. 2. We select the latter best-fit point as our primary 'benchmark point' in parameter space from hereon.  Table 2. Charges for all the SM fields at our benchmark point in parameter space (a e , a µ , a τ , a Y ) = (0.59, 1, 0, 0.87), which we obtained by fitting for (a e , a Y ) subject to a µ = 1, a τ = 0, and using all observables in the fit. This corresponds to the bottom line in Table 1. elliptical region in parameter space is selected. From hereon, we shall choose this best-fit point as our primary 'benchmark' for further study, at which we shall consider the impact of other important experimental constraints in §5. The best-fit values for a e , a Y , and the overall normalisation C are recorded in Table 1.
We can draw some interesting conclusions from these global fits, particularly from the fit to all observables which appears to select a clear preferred region. Firstly, while we see that a reasonable fit can be obtained with new physics only in the muon (which corresponds to the point a e = a Y = 0, which lies within the 2σ contour), there is strong pull to include some new physics component in the electron. Furthermore, independently of a e , we see that the fit also favours turning on a significant component a Y > 0, which we know gives a flavour-universal contribution to C 9 and C 10 .

Other phenomenological constraints
In addition to fitting the NCBAs, there are several important phenomenological constraints on our model. These come from high-p T LHC searches for the Z (e.g. in pp → µµ), B s meson mixing, a number of electroweak precision observables (in particular, LFU precision measurements in the Z boson couplings from LEP, and the measurement of the ρ-parameter), and neutrino trident production. In this Section we shall consider each of these constraints in turn.
In the same spirit as above, our goal here is not to characterise the phenomenology of these models in all detail, but rather to analyse the general structure of the different constraints within our framework, and to show as a 'proof of principle' that there is a viable parameter space. Thus, we shall present the theoretical expressions relevant to each constraint and comment on their form, and we are content to extract the numerical bounds only at our benchmark point in the parameter space. Recall that this benchmark point corresponds to the second line of Table 1, for which the charges are listed in Table 2.

Direct searches at the LHC
The strongest Z bounds from direct searches at the LHC can be obtained by interpreting recent ATLAS searches for resonances in both the pp → µ + µ − and pp → e + e − channels in 139 fb −1 of 13 TeV data [71], which extend up to a dilepton invariant mass of 6 TeV. 18 The constraints on a number of simple Z models coming from these ATLAS searches have recently been computed in Ref. [75], and if we wished to calculate an accurate bound from direct searches valid for any value of the coupling g X , then we should follow a similar methodology.
For our purposes in this paper, we first note that if m Z ≈ g X v Φ exceeds 6 TeV, then the constraints from these direct bump searches do not apply directly. In this region of parameter space, which we shall refer to as the 'contact régime', the high-p T tail of the dimuon invariant mass-squared distribution is of course still sensitive to the Z , but its effect may be computed using a simple EFT calculation involving the four-fermion effective operators which couple the final state lepton pair to a quark pair in the initial state [76]. We must first establish that this régime is even accessible within our framework; in other words, we must show that the Z mass can be as heavy as 6 TeV.
As we shall see in §5.2, there is an upper bound on the vev v Φ from B s mixing of approximately v Φ 4 TeV; there is also an upper bound on the coupling g X , beyond which the Z becomes strongly coupled and perturbativity breaks down, which we shall now estimate. This bound is approached when the width Γ Z of the Z resonance becomes broad, which we take to be when Γ Z /m Z ≈ 0.3 or so. To calculate Γ Z , we need the partial width of the Z decaying into a pair of massless fermions ff , which in the limit where m Z 2m t is given by Γ ff Z = Cg 2 XQ 2 f m Z /(24π), where C = 3 for quarks and C = 1 for leptons, andQ f denotes the U (1) X charge of the fermion f , as recorded in (2.4), and, at our benchmark point in parameter space, in Table 2. Summing over all the SM fermion species, 19 we obtain and so perturbativity breaks down for couplings as large as g X ∼ 0.3/0.071 ≈ 2.1. We may therefore consider large-ish couplings in the window 1.5 g X 2.1 for which the contact 18 There are also less constraining Z searches from ATLAS on smaller data sets in other channels, such as a pair of Z → tt searches in 36.1 fb −1 of 13 TeV data [72,73] which extend out to 5 TeV, and a Z → τ + τ − search in 10 fb −1 of 8 TeV data [74] which extends out to 2.5 TeV. 19 Other decay modes of the Z , for example to ZH, can be shown to be subleading.
régime is accessible, with the Z couplings remaining perturbative. In this régime we can compute the bound on v Φ using an EFT approach. This is in some sense the weakest possible bound on v Φ coming from direct searches, at least at our benchmark point in parameter space. In Ref. [76], bounds are tabulated for each semi-leptonic four-fermion operator using such an EFT approach, considered one operator at a time, which are valid in this contact régime. 20 For the benchmark point we are considering however, there are of course multiple relevant four-fermion operators turned on, with the dominant couplings due to the following three operators (we adopt the normalisation of Wilson coefficients used in Ref. [76] for ease of comparison): 2) Within our framework, the Wilson coefficients are given by Using the most recent 139 fb −1 ATLAS search described above [71], and including the contributions to the pp → + − cross-section from all four-fermion operators, one obtains the bound v Φ > 3.1 TeV (5.4) at the 95% C.L. 21 We shall see that this provides one of the most important bounds on the model at the benchmark point, alongside bounds from the ρ-parameter and from B s mixing which we compute next.

Neutral meson mixing
The quark flavour violation that mediates b → s transitions, which is a necessary ingredient for explaining the NCBAs, immediately results in a BSM contribution to B s meson mixing. Within our framework of Z models there is a contribution from tree level exchange of the Z , in addition to a contribution from scalar box diagrams involving Φ. This is exactly analogous to the meson mixing constraints derived in Ref. [31], which we here follow. 22 The mixing amplitude M 12 for the B s meson system takes the form Note that the bounds in Ref. [76] were computed using an ATLAS search on 36 fb −1 of 13 TeV data from 2017 [77], which was superseded by the 2019 search published in Ref. [71]. 21 We are very grateful to David Marzocca for performing this calculation and providing us with the result. 22 Given our simplifying assumption Y Qd = 0, there are no analogous bounds from kaon or B d mixing. There are nonetheless bounds from D mixing, though these are less constraining than those from Bs mixing, since the coefficient of the latter is inextricably tied to fitting the NCBAs. In particular, we find that D mixing constraints are much weaker than the ones from Bs mixing as long as there is a slight hierarchy in the couplings, viz. |YQs| < |Y Qb |.
where m W is the mass of the W boson and S 0 2.3 is a SM loop function. The Wilson coefficient C bs LL is given by (5.6) Here, the first term on the right-hand-side is due to the Z exchange, which scales somewhat unusually like v 2 Φ ; the two powers of v Φ in the denominator from the Z propagator are compensated by four powers of v Φ in the numerator arising from the square of the coupling The second term on the right-hand-side is due to the 1-loop box diagram.
While the mass difference ∆M s ∝ |M 12 | is measured with excellent precision, ∆M exp s = (17.757 ± 0.021)/ps [78], the SM prediction comes with a sizable uncertainty. Recent work using sum rule calculations of hadronic matrix elements quotes ∆M SM s = (18.5 +1.2 −1.5 )/ps [79]. If we instead use the lattice average of hadronic matrix elements from [80] (see also [81,82]) as well as |V cb | = (39.9 ± 1.4) × 10 −323 we find ∆M SM Within our framework, we see that the reasonable agreement of the SM prediction for B s mixing with the data, which precludes too large a BSM contribution to M 12 , provides an upper bound on the parameter v Φ , which cannot therefore be pushed arbitrarily high. 24 Given the constraints from direct searches place a lower bound of 3.1 TeV on v Φ (see §5.1), we have bounds squeezing the parameter v Φ from both sides. For the purposes of this paper, we seek only to show that there is a viable window of parameter space for v Φ between these bounds, for which it is sufficient to consider our benchmark point in parameter space (see Table 2). The B s mixing constraints on v Φ at the benchmark point are shown in Fig. 3, in the plane of Re(Y Qb Y * Qs ) vs. m Q . In these plots, the shaded green regions show the 1σ and 2σ best fit regions to the NCBA data, using all the observables, as discussed in §4. Note that the overall normalisation of the Wilson coefficients, which is extracted from the fit to the NCBAs, is proportional to Y Qb Y * Qs /m 2 Q from (4.4, 4.5). Hence the green regions correspond to bands with approximately fixed gradient (given the log scale of the plot). Each of the dashed contours then show the upper bounds on v Φ coming from B s mixing. Thus, if we take the central value of the fit to the NCBAS using all observables, we see that the maximum value of v Φ is between 3 and 3.5 TeV, depending of the precise B s mixing bound that is 23 This value is a conservative combination of the inclusive determination of |V cb | quoted in [1] and the two recent exclusive determinations from [83] and [84]. The sizable discrepancy between the inclusive and exclusive values is taken into account by rescaling the uncertainty by a factor 2.6, following the PDG prescription. 24 Note also that the right-hand-side of (5.6) does not depend on any of the parameters from our U (1)X gauge sector (i.e. gX , ae, aµ, aτ , or aY ), but only on the parameters from the mixing and scalar sectors.  imposed. If we fit the NCBAs at the 1σ (2σ) contours, we can loosen these constraints to v Φ 4 TeV (v Φ 5 TeV) respectively.
In the dark grey region in the upper left corner, B s mixing is never in agreement with measurements for any value of v Φ , because the (v Φ -independent) 1-loop contribution to (5.6) saturates the bound on its own. Thus, at the benchmark point the model is tightly squeezed by the combination of constraints on v Φ from direct searches and B s mixing, but there remains a viable window of unexcluded parameter space. Of course, one expects that by deviating from the benchmark point this viable window can be widened (or narrowed); however, a comprehensive analysis of the parameter space is beyond the scope of this paper.

Electroweak precision observables
Here we discuss two important constraints coming from electroweak precision observables (EWPOs), namely measurements of LFU in the Z couplings, and a constraint from the ρparameter. In general, these constraints from EWPOs arise due to the coupling of the Higgs to U (1) X , and so can be eased by dialling down the value of the parameter a Y which sets thê Q H ; while a Y is not essential for fitting the NCBA data, it is, interestingly, the parameter that determines the size of the axial component (i.e. the Wilson coefficient C α 10 ) of the flavour anomalies, which recall is necessarily flavour universal within our framework.

LFU of Z couplings
As we discussed in §3.2, the Z contains a small admixture of the U (1) X gauge boson X provided a Y = 0, and thus inherits some flavour non-universality in its couplings. Flavour non-universality in the leptonic decays of the Z are constrained by the LEP measurement [1] In the models we are considering, the ratio of partial widths is where g f f Z is the coupling of the physical Z boson to the fermion pair ff . One can obtain the couplings g f f Z by first writing down the terms in the lagrangian which couple the charged leptons to the neutral bosons B, W 3 , and X: We then rotate to the mass basis, and from X = sin α z Z +cos α z Z , this results in Z couplings that are suppressed by sin α z . To leading order in sin α z , we find the couplings are: where κ L ≡ − 1 2 g cos θ w + 1 2 g sin θ w and κ R ≡ g sin θ w (5.12) correspond to the couplings in the SM. We expand R model to leading order in sin α z to find After substituting in (3.16, 3.17) for sin α z and m Z , and the central experimental values g = 0.64 and g = 0.34, this becomes: Comparison with the LEP limits at the 95% C.L. yields the bounds: Of course, the bound vanishes when either a e = a µ (in which case the Z , and thus the Z also, couples universally to electrons and muons), or when a Y = 0 (in which case the Higgs becomes uncharged, and so there is no Z − Z mixing). At our benchmark point (see Table 2) we have (a e , a µ , a τ , a Y ) = (0.59, 1, 0, 0.87), and thus (a µ − a e )a Y = 0.357 > 0, so comparing with the first bound in (5.15) yields the constraint v Φ > 1.03 TeV. (5.16) This bound is therefore much weaker than that from direct searches computed in §5.1.

The ρ-parameter
The Z − Z mixing, which occurs at points in our parameter space where a Y = 0, also alters the mass of the Z boson away from the SM prediction. In particular, the ρ-parameter is no longer equal to unity at tree-level when a Y = 0. Global fits to electroweak precision data give the experimental constraint [1] ρ = 1.00039 ± 0.00019. (5.17) The diagonalisation of the neutral gauge boson mass matrix carried out in §3.2 gives a smaller eigenvalue m Z , 25 which is identified with the Z mass, that satisfies [85] where recall the mixing angle is sin α z = a Y g X √ g 2 +g 2 (m Z /m Z ) 2 + O((m Z /m Z ) 4 ), and the W boson mass is m W = vg/2 as in the SM. This results in the constraint 4.4 TeV (5.19) at the 95% C.L. limit. At the benchmark point in parameter space (a Y = 0.87) this constraint is in fact very aggressive, implying v Φ > 3.9 TeV, (5.20) even stronger than the constraint from direct searches -though still consistent with the 1σ best fit to the NCBAs (see Fig. 3). Of course, this constraint becomes weaker as one deviates away from the benchmark point we have studied in the direction of decreasing a Y .

Neutrino trident production
The neutrino induced production of a di-lepton pair in the Coulomb field of a heavy nucleus, a.k.a. neutrino trident production, is known to be an important constraint on models with gauged muon number [86]. Muonic tridents, ν µ → ν µ µµ have been measured at the CCFR experiment [87]. The quoted experimental value for the cross section corresponds to a 95% C.L. limit of The modifications of the trident cross section at the CCFR experiment coming from the exchange of a virtual Z gauge boson can be written as [88] σ where the new physics corrections to the effective couplings ∆g V and ∆g A are in our model given by At our benchmark point, the CCFR measurement given above implies the bound v Φ > 0.27 TeV. (5.25) This bound is considerably weaker than the ones from electroweak precision observables and from direct searches at the LHC.

Conclusions
We have proposed a new framework of Z models based on gauging an almost arbitrary linear combination of the accidental U (1) symmetries of the SM, i.e. baryon number and individual lepton numbers, as well as global hypercharge. Within this framework of models, the new physics associated with the Z at the TeV scale respects these global symmetries of the SM, whose breaking is therefore postponed until some even higher energy scale Λ LFV (at which neutrino masses are generated). Such a scenario is hinted at by the recent observations of possible new physics in rare B meson decays, since all these anomalous measurements respect the accidental symmetries of the SM. The conservation of lepton flavour, in particular, is naturally linked in these models to a violation of lepton flavour universality at the scale Λ LFUV at which the Z resides (since non-universal lepton charges are required to align the weak and mass eigenbases for charged leptons, which protects lepton flavour). Such a Z model therefore offers a tempting explanation of the neutral current B anomaly data, in which lepton flavour universality is observed to be violated between e and µ in b → s transitions, and we explore this possibility. The requisite quark flavour violation is introduced through a heavy vector-like quark state.
In the rest of the paper we explored the phenomenology of these models. Our analysis of the parameter space is not comprehensive. Rather, we were content to point out some interesting features that are common to all these models, such as the freedom to add a flavour-universal axial component to the B anomaly fit. We then made a number of wellmotivated assumptions to restrict the parameter space to a region of particular relevance to explaining the neutral current B anomalies, and by performing a global fit in this region we extracted a benchmark point in our space of models, at which to examine the phenomenology more carefully. At this benchmark point the Z couples to both left-handed and right-handed electrons and muons. We have shown that the model is consistent at this benchmark point with bounds from direct LHC searches, B s mixing, electroweak precision observables, and neutrino trident production.
Finally, in an Appendix we discuss how neutrino masses can be generated in such a setup, which involves a rather detailed construction of a dark sector. One can thence relate neutrino masses to the high energy scale Λ LFV at which the accidental symmetries of the SM are eventually broken within our framework. This gives an estimate Λ LFV ∼ 10 5 TeV, which is indeed much heavier than the scale Λ LFUV ∼ 1 TeV of new physics associated with the Z , and so is therefore consistent with our original hypothesis that Λ LFV Λ LFUV . An important next step would be to carry out a more comprehensive study of the phenomenological constraints on the parameter space of these models; for example, one might like to incorporate constraints from electroweak precision observables such as the ρ-parameter into the global fit to the flavour anomaly data, since we have seen that such precision observables do provide important constraints on our parameter space, in particular on the value of the parameter a Y .
Another promising future direction is to explore an alternative 'benchmark scenario' to the one we considered in this paper, in which the couplings of the Z to light quarks are set to zero, thereby significantly loosening up the constraints from direct searches at the LHC. This requires a Y = 0 and a e +a µ +a τ = 0, and so leaves two independent parameters to scan over in our fit to the flavour anomaly data. Given a τ is essentially unconstrained phenomenologically, one would most likely want to freely scan over values of a e and a µ (setting a τ = −a e − a µ ). Setting a Y = 0 also has the consequence of removing the tree-level Z −Z mixing, and so would also relax the bounds from electroweak precision observables (such as that coming from the ρ-parameter), which we found to be stringent constraints for large values of a Y . We therefore expect the phenomenology to be much more open in this region than in the benchmark we chose to study in this paper. Of course, taking this limit has significant implications for the flavour anomalies, since it also removes the possibility of any axial component in the flavour anomaly explanation. This scenario gives a simple illustration of the interesting interplay between different experimental constraints and the character of the flavour anomalies within our framework of gauging the accidental symmetries of the SM. bounds on our model (at the benchmark point in parameter space) coming from pp → µ + µ − and pp → e + e − direct searches, described in §5.1. We also thank other members of the Cambridge Pheno Working Group for helpful comments. JD is supported by The Cambridge Trust and by the STFC consolidated grant ST/P000681/1. The research of WA is supported by the National Science Foundation under Grant No. PHY-1912719. WA acknowledges support by the Munich Institute for Astro-and Particle Physics (MIAPP) of the DFG cluster of excellence "Origin and Structure of the Universe".

A Neutrino masses and the dark sector
Recall that a primary goal of this paper was to extend the SM in a way that preserves its global symmetries, in particular U (1) B and each individual lepton number U (1) L i , which are experimentally tested to very high precision. From this premise, we arrived at the fourparameter family of anomaly-free U (1) X charge assignments recorded in (2.4), for which the alignment of the charged lepton mass basis with the weak eigenbasis leads to a natural protection of lepton numbers, while at the same time predicting lepton flavour universality violation.
Nonetheless, we know that the lepton number symmetries cannot, ultimately, be exact symmetries of Nature, due to one important observation: the measurement of neutrino oscillations. Within the framework we have put forward in this paper, there is (by design) no lepton flavour violation in charged leptons, at least at the renormalisable level to which we have worked so far, and so we are led to suppose that all the observed lepton flavour violation (as parametrized, for example, by the PMNS matrix elements) must originate from the neutrino sector. Recall that in addition to the left-handed weakly-interacting neutrinos of the SM there are three right-handed SM singlet states in our setup, which we introduced to 'soak up' the U (1) 3 X and gravitational anomalies; the natural interpretation of these dark states is to identify them as right-handed neutrinos, which are charged only under the U (1) X gauge symmetry.
The measurement of neutrino mass-squared differences of the order 10 −3 eV 2 does not, however, point us ambiguously towards an energy scale Λ LFV (to use the terminology introduced in §1) at which the lepton number symmetries become broken, because such a cut-off scale is highly dependent on the physical mechanism which gives rise to the neutrino mass terms. However, with the charge assignments presented in (2.4), we in fact find a neutrino mass sector that is in tension with our hypothesis that Λ LFV resides much higher than the TeV scale (Λ LFUV ) associated with the Z boson and corresponding LFUV effects, where Λ LFUV ∼ v Φ . In this Section we shall explain why there is such a tension in the neutrino sector, before sketching how the dark sector of the model can be altered to resolve this tension, and thus give a compelling account of neutrino masses and the eventual breaking of lepton number symmetries.
The reason for the apparent tension is as follows. Firstly, one can use the three dark states ν i R to write down a diagonal matrix of renormalisable Yukawa couplings for neutrinos, of the form 26 We can also write down Majorana mass terms involving the SM singlet states ν i R , possibly with insertions of the scalar field Φ to soak up the U (1) X charge, depending on the particular values of a e , a µ , and a τ . For example, in the benchmark case that we defined in §4.2 (and have studied phenomenologically in §5), we can write down the following operators where M indicates a mass scale which is a priori unrelated to any other scale in the theory, and the couplings α, β, and γ are all dimensionless coefficients. After spontaneous breaking of U (1) X , the dimension four operators in (A.2) lead to Majorana mass terms with coefficients set by the scale v Φ , which we know is of order a few TeV. These effective Majorana operators break lepton flavour symmetries, in this case U (1) µ and U (1) τ , at the scale of v Φ . Now, if this lepton number violation were confined to the dark sector states this would not pose the problem, since from the point of view of the SM fields the SM's accidental global symmetries would remain intact up to the higher scale Λ LFV . But that is not the case, because the dark states ν i R necessarily interact with the charged leptons i through the diagonal Yukawa interactions at a low energy scale, 27 and we therefore have no natural mechanism for preserving lepton flavour symmetries in the SM sector up to the high scale Λ LFV . Thus, the tension arises precisely because of the interplay between the lepton flavour violating Majorana interactions with the Yukawa couplings which couple the dark sector to the SM.

A.1 Alternative dark sectors
It is possible to resolve this tension by exploring alternative dark sectors. 28 To see how this can be done, recall that the dark sector states ν i R were introduced to 'soak up' the nonvanishing gauge (and gauge-gravity) anomalies. We introduced three dark states because that was sufficient to guarantee cancellation of all anomalies for any choice of the four rational coefficients (a e , a µ , a τ , a Y ) -see equations (2.6, 2.7).
We here remark in passing that whatever extra chiral states we introduce the 'soak up' any field theory anomalies are better off being dark, for two reasons. Even if we added SM non-singlet states that were chiral only under U (1) X , these states would be on the one hand dangerous from the point of view of phenomenology, since they would acquire masses set by 26 Of course, even if the Yukawa couplings yi were tiny, these terms on their own cannot explain the neutrino mass sector, because the PMNS matrix features large mixing angles. 27 One rather unsatisfying resolution to this problem is that there is in fact a discrete Z2 symmetry which bans these Yukawa interactions. 28 We note in passing that the idea of using extra chiral states which 'soak up' gauge anomalies to serve as dark matter has been explored before in the context of gauging U (1)B−L [89][90][91][92]. In this paper, we do not consider dark matter phenomenology. v Φ ∼ 3 TeV multiplied by Yukawa-like dimensionless couplings (and so might easily have masses of 1 TeV or lower), for which the experimental bounds on, say, coloured states are strong (see for example recent bounds from gluino searches [93]), and on the other hand less useful from the point of view of soaking up anomalies, since (for example) an SU (3) triplet would always appear in (2.6) with a multiplcity factor of 3. If the added states were completely chiral under the SM the situation is worse still, because their masses would then be set by the Higgs vev v, which is very problematic for LHC phenomenology even if such states were colourless. However, as we noted above (see footnote 12), the assignment of dark charges in (2.7) does not provide the only solution to the Diophantine equations (2.6). As we there observed, for particular values of (a e , a µ , a τ ) there may be other solutions for the dark charges which soak up the anomalies; generically, these other solutions will not allow Yukawa couplings of the form (A.1), and in this way may circumvent the tension presented above. We will at times refer to these other solutions as 'non-trivial solutions' to (2.6), since the choice made in (2.7) is trivial from the number theory perspective. Indeed there are even non-trivial solutions with only two dark states, 29 and certainly there will be many more if we allow say four or more dark states. The general solutions presented in Ref. [69] can be adapted to generate solutions to this pair of equations at will, for any number (greater than one) of chiral dark states.
To furnish us with a concrete example of such an alternative dark sector, consider the benchmark case we set up in §4.2. By normalising to a µ = 1 and setting a τ = 0 for simplicity, 30 we found the best-fit point by performing a global fit to the NCBA data (see Fig. 2), where we assumed the anomalies in (2.6) were soaked up by states with chargesQ ν 1 = 0.59,Q ν 2 = 1, andQ ν 3 = 0 (i.e. exploiting the 'trivial solution' to the anomaly equations). But now consider an anomaly-free model with an alternative dark sector, also with three dark chiral fermion states, in which (a e , a µ , a τ ,Q ν 1 ,Q ν 2 ,Q ν 3 ) = corresponding to a non-trivial rational solution to the anomaly cancellation equations (2.6). This alternative anomaly-free charge assignment coincides (very nearly) with the benchmark case when restricted to the SM fields, and so shares the same phenomenology. The dark sector is however very different, which leads to a different story concerning neutrino masses, which we shall soon explain. 29 For example, with only two dark states ν 1 R and ν 2 R , there are anomaly-free solutions with (ae, aµ, aτ , −Q ν 1 , −Q ν 2 ) equal to (any permutation of) the sets (7, 8, −9, −1, −5), (9, −2, −10, −4, 7), etc. 30 We emphasize that the choice aτ = 0 was really just to simplify the discussion and the notation, rather than simplify the physics; the value of aτ does not affect any of the relevant phenomenological bounds computed in §5.
The crucial point is that with these dark sector charges we can no longer write down renormalisable Yukawa couplings of the form (A.1) which couple the dark sector to the SM fermions. At the renormalisable level, the neutrinos are now strictly massless, and we must pass to the SMEFT to explain the origin of neutrino masses.
Before we do so, we should first give masses to the dark states. While massless dark fermions of this kind are not in conflict with data from collider experiments, they would have profound consequences for the cosmological evolution of the Universe. To avoid a detailed analysis of the cosmological constraints, one can make the dark fermions heavy by introducing a pair of extra dark scalars transforming in representations where here χ a indicates either of the dark scalars in (A.5). M dark is a rank-3 mass matrix and so leads to three non-zero masses for the dark fermions, all at the scale of U (1) X breaking (i.e. the TeV scale).

A.2 Neutrino masses: estimating the scale of lepton flavour violation
If for simplicity we take a Y = 9/10 (which is close to the best-fit value of a Y = 0.87 obtained in §4.2), the left-handed lepton doublets have the rational chargeŝ which are, by construction, numerically almost equal to the charges in the benchmark case as recorded in Table 2.
One can show that mass terms involving left-handed neutrinos first appear at dimension six in the SMEFT, due to lepton flavour violating operators of the form where again χ a here denotes any dark scalar charged under U (1) X . After spontaneous breaking of U (1) X , these dimension six operators reduce to the familiar dimension five Weinberg operators. Note the appearance of the scale Λ LFV in this EFT expansion; recall Λ LFV reflects the scale of new physics which may break the global symmetries of the SM, in particular lepton flavour, in other words the scale at which our Z model breaks down. The lower scale Λ LFUV ∼ 1 TeV is the cut-off scale at which the SMEFT breaks down, being resolved at short-distances by our Z model. We can populate the mass matrix c ij up to a single zero, viz.
if we introduce three more dark scalars transforming in representations which is known to be a sufficiently dense texture to accommodate present data on neutrino masses and mixings. We emphasize that our construction of a viable dark sector in this Appendix is intended as a proof of principle, and we do not intend for this setup to be understood as being in any way a 'minimal choice'. 31 Finally, we can now estimate the higher cut-off scale Λ LFV associated with lepton flavour violation, using the scale of neutrino masses. At the level of naïve dimensional analysis we require v 2 v Φ Λ 2 LFV ∼ m ν ≈ 10 −13 TeV =⇒ Λ LFV 10 5 TeV, (A.11) which sure enough far exceeds the TeV scale associated with the Z and LFUV effects such as the measurements of R K ( * ) .
B Details of the fit to rare B decay data In this Appendix we provide details about the fit to rare B meson decay data that we perform to identify the preferred parameter space of our model. We carry out two fits in the 3-dimensional parameter space of Wilson coefficients C µ 9 , C e 9 and C µ 10 = C e 10 ≡ C 10 using flavio [67]. In fit 1, we include the measurements of LFU ratios R K and R K * from LHCb [3,4] and Belle [94,95], the combination of the B s → µµ branching ratio from [13], that includes data from LHCb, CMS, and ATLAS [5][6][7][8], as well as measurements of the branching ratios of the inclusive decays B → X s from BaBar [96] and Belle [97]. We find a compact region of parameter space that is compatible with the considered data and we approximate the best fit region by a multivariate Gaussian. The corresponding central values for the Wilson coefficients, their uncertainties and the correlation matrix are In fit 2, we include in addition measurements of the branching ratios of the exclusive semileptonic decays B ± → K ± µµ, B 0 → K 0 µµ, B 0 → K * 0 µµ, B ± → K * ± µµ, B s → φµµ and Λ b → Λµµ, as well as all available measurements of CP averaged angular coefficients in these decays from LHCb, CMS, and ATLAS [10][11][12][98][99][100][101][102][103][104]. The best fit region is given by  .7), the best fit regions in the parameter space of the Wilson coefficients can be mapped onto the model parameters a e , a Y , and C. The result is shown in Fig. 2 for the a e vs. a Y plane, profiling over the normalisation factor C.