Complete Vector-like Fourth Family with U(1)′: A Global Analysis

In this paper we present an in-depth analysis of a recently proposed Standard Model extension with a complete fourth generation of quarks and leptons, which are vector-like with respect to the Standard Model gauge group and charged under a new spontaneously broken vector-like U(1)′ gauge symmetry. The model is designed to explain the known muon anomalies, i.e. the observed deviations from Standard Model predictions in the anomalous magnetic moment of the muon, ∆aμ, and in b → s`+`− processes. We perform a global χ2 analysis of the data with 65 model parameters and including 98 observables. We find many points with χ2 per degree of freedom ≤ 1. The vector-like leptons and the new heavy Z ′ are typically much lighter than a TeV and would, thus, be eminently visible at the HL-LHC. ∗kawamura.14@osu.edu †raby.1@osu.edu ‡trautner@mpi-hd.mpg.de ar X iv :1 91 1. 11 07 5v 1 [ he pph ] 2 5 N ov 2 01 9


Introduction
In the search for physics beyond the Standard Model (SM) one certainly looks for new particles produced at the LHC. Absent any new discoveries, one then considers experimental discrepancies with SM predictions. There will always be 2σ or 3σ discrepancies, some of which are real and some of which are just statistical fluctuations. In an attempt to distinguish these two possibilities, one might focus on a cluster of discrepancies which can all be resolved with the same new physics. This is what we have done in a previous letter [1]. We have shown that the muon anomalies associated with the anomalous magnetic moment of the muon, ∆a µ , and the angular and lepton non-universality anomalies in b → s + − decays can simultaneously be resolved by the addition to the SM of a complete vector-like (VL) family of quarks and leptons together with a new U(1) gauge symmetry carried only by the new VL states. VL leptons are introduced for ∆a µ in Refs. [2][3][4][5]. The b → s + − anomaly is addressed by introducing a Z boson [6][7][8][9] and new particles which induce box diagram contributions [10][11][12][13][14][15][16][17][18]. Both anomalies can be explained simultaneously in models with VL fermions and a Z boson [19][20][21][22][23].
In the present paper, we perform a global χ 2 analysis of these phenomena with the addition of all SM processes that might feel the effects of mixing between the VL family and the SM families. We find many solutions with χ 2 /N d.o.f. ≤ 1. Moreover, the model is highly testable via both direct observation of new physics at the LHC or via the improved analysis of SM processes. For example, precision measurements of K-K and B d -B d mixing may be sensitive to the new physics. Finally, the tight observational constraints on the µ → eγ branching ratio hinders a simultaneous fit to ∆a µ and ∆a e . We can only fit one but not both, and we choose to fit the former.
The paper is organized as follows. In Section 2 we briefly review the details of the model and describe the mass mixing between the SM and VL states. In Section 3 we provide the theoretical formulae for calculating the 98 observables in our analysis. In Section 4 we present the results of the χ 2 analysis with many plots illustrating the range of VL masses and the quality of the fits. In addition we choose four "best fit points" to illustrate some of the new physics processes which are, in principle, observable at the LHC. Finally, in Section 5 we summarize our results. More details of the fits are presented in the Appendices.

Matter Content and Fermion Mass Matrices
We study a model with a complete VL fourth family and U(1) gauge symmetry. The quantum number of all particles are listed in Tables 1 and 2 3)

4)
(2.9) For electrically charged fermions, f = u, d, e, the mass basis is defined aŝ Here, m Ea , m Ua , m Da (a = 1, 2) are masses for the new charged leptons, up and down quarks, respectively. These are predominantly the VL fermions of the gauge basis. We consider the type-I see-saw mechanism to explain the tiny neutrino masses. Since the U(1) gauge symmetry prohibits Majorana masses for the VL family, only three families of right-handed neutrinos have Majorana masses, (2.14) The Majorana masses are assumed to be M ij Maj ∼ 10 14 GeV. The mass matrix for the neutrinos is then a 10 × 10 matrix, where the Dirac mass matrix M n is defined in Eq. (2.7). The 5 × 5 Majorana mass matrix is given by (2. 16) After diagonalizing this matrix, there are three left-handed Majorana neutrinosν L i with mass of O (v 2 H /M Maj ), and three right-handed Majorana neutrinosν R i with mass of O (M Maj ) as in the ordinary type-I see-saw mechanism. In addition to these, there are two Dirac neutrinos N 1 , N 2 with mass of O (v φ ) which are predominantly the VL neutrinos of the gauge basis. Mixing between the left-and right-handed neutrinos is suppressed by the huge Majorana mass terms. An approximate mass basis is then defined as 2 n L := (U n L ) † n L ,n R := (U n R ) † n R , (2.17) with unitary matrices U n L and U n R , given in Appendix B where we discuss the diagonalization of the neutrino mass matrix more explicitly.
There are three electrically neutral scalar fields in this model. We expand them around their VEVs as Here, h, χ, and σ are physical real scalar fields while the pseudo-scalar components a H and a χ are eaten by the gauge bosons. We introduce effective quartic couplings of the scalars χ and σ, which parametrize the scalar masses as (2. 19) In this paper, we do not specify a scalar potential in this model and treat their VEVs and the effective quartic couplings as input parameters. For an effective quartic coupling λ χ in a perturbative regime and a sizable gauge coupling g , the masses of χ and the Z gauge boson should roughly be of the same order. This is important, since the Yukawa couplings of χ together with v Φ set the scale of mass mixing between VL particles and the SM families. The fact that this mixing is necessary for an explanation of the muon anomalies means that also χ will give sizable contributions in our fits. In contrast, v φ could be very large, and therefore m σ very heavy, as long as the φ-Yukawa couplings are small enough to prevent a complete decoupling of the VL fermions. The scalar σ, thus, can be very heavy and its contributions irrelevant for all discussed observables. In fact, this limit resembles the case of rigid tree-level VL masses. Consistent with that, contributions from σ are always negligible at the best fit points we discuss below.

Gauge and Yukawa Couplings
In the gauge basis, there are no interactions between the Z boson and the SM fermions. In order to explain the muon anomalies, the SM families are required to mix with the VL families. Consequently, also the SM gauge couplings of quarks and leptons will receive corrections from these mixing effects. Of course, these couplings must remain consistent with all SM observables and we shall verify this. For this discussion we combine left-and right-handed fermions to Dirac fermions as (2.20)

Neutral Gauge Couplings
The fermion couplings to the Z boson are given by where the charge matrices in the gauge basis are The coupling matrices in the mass basis arê The Z boson couplings are given by where Q f and T f 3 are the electromagnetic charge and third component of the weak isospin for a fermion f , respectively, while s W (c W ) denotes the (co)sine of the weak mixing angle. The flavor space projectors are defined as P 5 := diag (0, 0, 0, 0, 1) and P 5 := 1 5 − P 5 = diag (1, 1, 1, 1, 0) .

(2.27)
The coupling matrices in the mass basis are given bŷ (2. 28) In general, this model has tree-level flavor changing neutral vector currents mediated by the Z boson. However, for the SM generations these are automatically suppressed by O(m 2 f SM /M 2 F VL ) coefficients which we show analytically in Appendix B.
Here, m f SM and M F VL denote the mass of the SM fermion f SM , as well as the mass of the VL fermion F VL .

CKM and PMNS Matrices
The fermion couplings to the W boson are given by W + µ [uγ µ (P 5 P L + P 5 P R ) d + nγ µ (P 5 P L + P 5 P R ) e] + h.c. (2.29) = W + µ ûγ µ ĝ W q L P L +ĝ W q R P R d +nγ µ ĝ W L P L +ĝ W R P R ê + h.c. (2.30) In the mass basis, the coupling matrices arê (2.32) The extended 5 × 5 CKM matrixV CKM is defined aŝ Since one of the left-handed quarks is an SU(2) L singlet this extended CKM matrix has only rank 4 and is, therefore, clearly non-unitary. Correspondingly, there exist righthanded charged current interactions which are completely absent in the SM. Also the 3 × 3 sub-matrix, which corresponds to the SM CKM matrix, is generally non-unitary due to mixing with the VL quarks. Again these effects are suppressed by O(m 2 f SM /M 2 F VL ) coefficients. In addition, the right-handed current interactions to the SM quarks are also suppressed and at most O (10 −4 ) as shown in Eq. (B.64), see Appendix B. These are negligible against experimental sensitivities.
The mixing between the SM and VL neutrinos are suppressed by the huge Majorana mass term. In Appendix B, we found that non-unitarity of the 3 × 3 PMNS matrix is induced by the Yukawa coupling λ n and tiny mixing angles between the SM and VL charged leptons, that is U 2 e L in Eq. (B.18). In other words, there would be non-zero mixing between the SM neutrinos and the Dirac neutrinos if λ n ∼ 1. This is an interesting possibility, but is beyond the scope of this paper. We consider a parameter space where λ n 1 and the PMNS matrix is approximately unitary. The details of neutrino mass diagonalization as well as the definition of the PMNS matrix are shown in Appendix B.

Yukawa Couplings
The Yukawa interactions with the real scalars are given by Here the Yukawa matrices for the fermions in the gauge basis are given by with the mass matrices of Eqs. (2.6)-(2.9). In the mass basis, these arê As in the Z and W boson couplings, the flavor violating couplings to the Higgs boson is

Landau Pole Constraints on the U(1) Gauge Coupling
As already discussed in more detail in [1], the U(1) gauge coupling g diverges at a Landau pole at the scale (2.38) Here µ Z ∼ 1 TeV is the scale where we define the couplings of our model. In order for the model to be consistent up to a typical GUT scale of Λ g ∼ 10 16 GeV, we require g (1 TeV) < 0.35 in our analysis.

Observables
In our model, ∆a µ is explained by chirally enhanced 1-loop corrections involving the Z boson and VL leptons. At the same time, tree-level Z exchange induces new contributions to b → s + − . An explanation for ∆a µ requires sizable Z couplings to muons, in agreement with those necessary to explain deviations in b → sµµ. The required mixing of the SM and VL fermions may, thus, induce new physics effects in various observables in both, lepton and quark sectors. We have already shown that this model can explain the muon anomalies without spoiling other observables in Ref. [1]. The purpose of the present paper is to completely map out the parameter space where muon anomalies are explained consistently with all other observables, and at the same time discuss the consequences for future experiments. In this section, we introduce the 98 observables included in our χ 2 analysis. In our analysis, 1σ uncertainties for observables which only have upper bounds are set so that 1.64σ (1.96σ) deviation gives a value of 90% (95%) C.L. limit.

Charged Leptons
Here we study the charged lepton masses, lepton decays, as well as lepton anomalous magnetic moments. Central values and uncertainties of the observables are listed in Table 3.
The 1-loop beyond the SM corrections involving Z , Z and W bosons to ∆a µ are given by [4,37] Here, m e(n) B is the mass of the B-th generation charged(neutral) lepton. The flavor index B = 1, . . . , 5 runs over all five families, while b = 4, 5 runs only over the VL family. The loop functions are defined as The scalar 1-loop beyond the SM contributions to ∆a µ are given by [4,37] where S = χ, σ. Here, y h e b := m 2 e b /m 2 h and y S e B := m 2 e B /m 2 S . The loop functions are defined as F S (y) = − y 3 − 6y 2 + 3y + 6y ln (y) + 2 6(1 − y) 4 , (3.12) Analogous formulae for ∆a e are obtained by formally replacing 2 → 1 and µ → e. We now discuss the leading contributions to ∆a µ analytically. This will be important to understand upper bounds on v Φ and the masses of VL leptons stated below. The new physics contribution to ∆a µ is dominated by chirally enhanced effects proportional to λ e v H , namely the Higgs-VEV induced mixing between left-and right-handed VL leptons. Contributions involving the SM bosons are very suppressed. The leading contribution can be estimated as where s 2µ L(R) are mixing angles between the muon and VL leptons and C LR is a function of the mass ratios that we will define shortly. The mixing angles are approximately given by (see Appendix B for details) Here, M E L and M E R the masses of the doublet-and singlet-like VL leptons which can be approximated as The dimensionless function C LR is defined in (B.46) in Appendix B. It can be approximated as Figure 1 shows contour plots of the functions C Z and C χ . C Z has a maximum value ∼ 0.272 at is always positive (negative) and |C Z | > |C χ | at most parts of the parameter space. Altogether, an upper bound on ∆a µ is given by The equality is saturated when λ e ∼ −1.0, s 2µ L ∼ s 2µ R ∼ 1 and C χ ∼ 0. Inserting the maximal value of C Z one finds Thus v Φ 1.7 TeV is required to explain ∆a µ . Moreover, is required to maximize C Z . We are interested in upper bounds on the lightest VL charged lepton. For a fixed lightest VL charged lepton mass, the function C Z is maximized if the heavier state has the same mass, i.e. M E L = M E R := M E . Then x L = x R := x, and (3.20) Figure 2 shows contours of ∆a max µ = 2.68 × 10 −9 in the (m Z , M E ) plane where C Z is replaced by Eq. (3.20) and the gauge coupling constant g is fixed. Different colors correspond to different values of g . ∆a µ = 2.68 × 10 −9 can be realized only inside the contours for a given g . We further restrict the VL masses by The last inequality comes from our requirement for perturbativity, λ E 2 ≤ 1. This condition is depicted by the straight lines in Fig. 2. Altogether, the upper bound on the VL lepton mass is about 1.4 TeV where m Z ∼ 500 GeV and g = 0.35. Note that m E 1 ∼ 1.4 TeV is realized only if all of the conditions are satisfied: (1) 35. Consequently, the upper bound is hardly ever saturated. The dashed lines in Fig. 2 show the same contour but the destructive χ contribution with m χ = 100 GeV is added to C Z . m χ = 100 GeV is chosen to minimize the C χ contribution. Including the scalar contribution, the upper bound on  The dots and crosses represent best fit points (which in particular fit ∆a µ in the 1σ allowed region). Details will be explained in detail in Section 4. Red (green) dots have χ 2 < 28 (33), while points which are excluded by other data are shown as crosses or pluses. m E is tightened to 1.3 TeV. Clearly the actual upper bound will be lower if some of the conditions (1)-(6) are not satisfied. The points shown in Fig. 2 are results of our fit. As anticipated, good fits are only obtained for points within the contours. The details of our analysis will be shown in the next section.
On a different note, a lighter scalar χ allows one to explain ∆a µ with smaller Z contribution (see Figure 1), especially when the VL leptons are heavy. y L,R 1 or y L,R 1 is favored to suppress the destructive contributions from C χ . In the phenomenologically viable parameter space, M E L,R 250 GeV and a lower bound on y L,R is For y L,R ∼ 0.0625 with m χ = 1 TeV, we have C χ ∼ −0.06 . On the other hand, y L,R ∼ 100, corresponding to C χ ∼ 0.004, is possible if m χ ∼ 100 GeV and M E L,R ∼ 1 TeV. Another important consequence of ∆a µ is that the Higgs Yukawa coupling λ e should Hence, ∆a µ is proportional to λ e . For this reason, the U(1) charge assignment in our model must not be universal for (L L , E R ), but must be flipped as in Table 2. Importantly, such a charge assignment is incompatible with SO(10) unification. However, it is still consistent with unification in the Pati-Salam gauge group, The dominant decay modes of the charged leptons are three-body decays via a W boson. The branching fraction for a lepton i to decay into a lighter lepton j is given by where m i and Γ i are the mass and decay width of the lepton i , respectively. The function F is given by Experimental values of the lepton masses and decay widths are listed in Table 4. For the muon decay rate, the tree-level branching fractions are multiplied by a QED correction factor η QED = 0.995802 [25]. This factor is less important for tau decay. Just as for the charged lepton masses, the charged lepton decay rates are measured more precisely than the accuracy of our numerical analysis. We assume 0.01% relative uncertainties for these observables, remarking that we could always fit them by increasing the numerical accuracy of our analysis. Branching fractions could deviate from their SM predictions if the mixing between the SM families and VL family affects the W couplings. The values obtained in our model are compared with the tree-level SM values, that are given by

3.1.3
i → j γ Lepton Flavor Violating (LFV) processes are severely constrained by experiments. We follow [38] to calculate one-loop corrections including general gauge and Yukawa interactions. The Lagrangian for general gauge and Yukawa interactions is given by (3.24) where k are external charged leptons, F internal fermions, V vector bosons and S scalars. The gauge couplings in our model are identified as where A, B = 1, . . . , 5, while b = 4, 5 runs only over the VL neutrinos. 3 The Yukawa couplings are given by The branching fraction is then given by where Here, y i , y i , and k j (with i = 1, . . . , 4 and j = 1, 2, 3) are loop functions defined in Ref. [38] while combinations of couplings are defined as with B = S, V, W . σ R is obtained from σ L by formally replacing ρ BF ↔ λ BF and ζ BF ↔ υ BF .
Just as for ∆a µ , the dominant contribution to i → j γ is again a chirally enhanced effect. To a good approximation, σ L is given by where the loop functions G Z,S are the same as in Eqs. (3.7) and (3.13). Using the results from Appendix B, analytic expressions for the branching fractions of µ → eγ and τ → µγ are given by with C LR given in Eq. (3.16). Here, e L,R ( τ L,R ) are the mixing angles between the leftand right-handed electrons (tauons) and the VL leptons, respectively. e L,R 10 −6 and τ L,R 10 −2 are required to suppress µ → eγ and τ → µγ, respectively. Once both of these processes are suppressed, τ → eγ is automatically suppressed as well.

3.1.4
i → j k l The neutral bosons also mediate LFV three-body decays, such as µ − → e − e + e − , τ − → µ − µ + µ − and so on. Effective interactions relevant for a decay − The branching fraction is given by [39,40] where masses of daughter leptons are neglected. Table 3) are given by The branching fraction is given by [41] (3.38) In this model, the Wilson coefficients are given by where X, Y = L, R. These LFV three body decays are dominated by Z boson exchange. Using the result in Appendix B, the Z contributions to µ − → e − e + e − and τ − → µ − µ + µ − are estimated as where s µ , e , and τ are the maximum values of s µ L,R , e L,R and τ L,R , respectively. BR(µ → 3e) is strongly suppressed by 6 e and will be much smaller than BR(µ → eγ). On the other hand, BR(τ → 3µ) scales as 2 τ , and therefore in the same way as BR(τ → µγ). BR(τ → 3µ) BR(τ → µγ) is expected because the former is proportional to an absolute sum of different chirality structures, while the latter is dominated by the left-right mixing effect. All other τ decays are suppressed by additional factors of e .

EW Bosons
The fermion couplings to the SM bosons, namely Higgs, W and Z bosons, are also modified by the mixing to VL fermions. This might affect their decays. For instance, LFV Higgs boson decays are predicted in models with VL leptons studied in Refs. [4,20,42]. All observables here are calculated at tree-level, except for h → γγ. All formulae that we use to compute two-body decays are summarized in Appendix A. Table 5 summarizes experimentally determined values of constants used in the EW boson observables. Experimental central values and uncertainties of the relevant observables are summarized in Table 6.

W Boson Decays
There is a right-handed charged current coupling to the W boson which is absent in the SM. Furthermore, the non-unitarity of the CKM and PMNS matrices can affect the W boson couplings. These will alter W boson decays. The branching fractions for W boson decays are given by where The function λ is defined in Eq. (A.3). SM predictions for these decays are calculated by replacing Here, experimental absolute values of the PMNS and CKM matrix elements are used. For the leptonic decays, radiative corrections and experimental uncertainties are small. We BR(Z → eτ ) 0.0 5.0×10 −6 Ref. [25] BR(Z → µτ ) 0.0 6.1×10 −6 Ref. [25] A e 0.1469 1 % SM use a 0.1% relative uncertainty for these decays. For the hadronic decay modes, QCD corrections may change the values by a factor proportional to α s /π ∼ 0.038. We use a relative uncertainty of 3.8% for the total hadronic branching fraction, while we use a 10% relative uncertainty for W → cs because experimental uncertainties here are still large.

Z Decays and Asymmetry Parameters
The Z boson couplings are, in general, changed by the mixing effects. This can affect the branching fractions and asymmetry parameters of Z decays. The Z boson couplings depend on the weak mixing angle θ W . For the lepton couplings, we use the effective angle s = 0.23154 including radiative corrections, while the tree-level value s W = 0.22343 is used for the quark couplings [25]. Using the effective angle is necessary to reproduce the observed values of A =e,µ,τ . The branching fractions for flavor conserving decays are given by for a fermion f i . Those for flavor violating decays are given by For BR(Z → had), the experimental value of a four-body decay BR(Z → bbbb) = 0.00036 [25] is added to the sum of two-body decays to quarks. The Z-pole asymmetry parameters are defined as where f i = e, µ, τ, s, c, b, and (axial-)vector couplings are obtained as SM predictions for these observables are obtained by formally replacing P 5 , P 5 → 1 3 , 0 3 in Eq. (2.28). The leading QED and QCD corrections to these decays are proportional to 3α e /4π ∼ 0.0019 and α s /π ∼ 0.038 for leptonic or hadronic decays, respectively. We, therefore, use a 0.19% (3.8%) relative uncertainty on the leptonic (hadronic) decays. The relative uncertainties for the asymmetry parameters are taken as 1% for A e , A τ , A c , A b , and 10% for A µ and A s , consistent with their experimental uncertainties. The uncertainties for flavor violating decays are determined from their experimental upper bounds. Let us illustrate how the SM boson couplings, in general, are very close to their SM values. We show this analytically in Appendix B. In general, one finds the mixing to be suppressed by Considering, for example, the left-handed Z − µτ coupling we find that it can be estimated as Corrections for lighter flavors are even more suppressed.

Higgs Decays
The Higgs boson couplings to SM fermions can, in general, depart from their SM predictions due to misalignment of the Yukawa couplings and mass matrices. We have studied the signal strengths for the measured decay modes to µµ, τ τ , bb, and γγ final states as well as the branching fractions for the unobserved decays ee, eµ, eτ , and µτ . Central values and uncertainties are set to their experimentally observed values. Decay widths for flavor conserving decays are given by and those for flavor violating decays are In addition to the tree-level decays, the VL families may significantly contribute to the loop-induced decay, h → γγ and h → gg. The decay width for h → γγ is given by [43] Here, f A runs over all the fermions in this model and A = 1, . . . , 5 is a flavor index. N f c is the number of color degrees of freedom and f ] AA is a diagonal Yukawa coupling constant of the Higgs boson to a fermion f A . The form factors are given by Similarly, the decay width of h → gg is given by [43] Γ where q A only runs over the quarks. Naively, one expects the size of VL fermion contributions to these one-loop decays to be suppressed by the squared ratio of SU(2) × U(1) breaking mass to VL mass, i.e. by a factor (λ This is exactly what we find. Using the result of Appendix B, contributions from VL fermions F L , F R are given by Here, λ f v H M F L,R and λ f λ f have been assumed. For the last equality in the first line, we use the series expansion A H 1/2 (τ ) ≈ 4/3 + 14/45 τ around τ ≈ 0. A possible cancellation between the two VL fermions gives an extra suppression. Altogether, we confirm that the VL fermions only give very small corrections to these decay rates. Especially VL quarks will be heavy, and their effects therefore particularly suppressed. Thus, there is no meaningful constraint arising from h → gg for our analysis, and also the Higgs boson production rate is unchanged with respect to the SM.
The signal strengths of the Higgs boson are defined as

Quarks
We study the SM quark masses, 9 absolute values of the CKM matrix elements and 3 CP phases α, β, γ. The Wilson coefficients relevant to the b → s + − processes are fitted to explain the anomalies. The new physics contributions will also affect neutral meson mixing, namely K-K , B d -B d , B s -B s and D-D mixing, (semi-)leptonic decays of B mesons and top quark decays. Central values and uncertainties for quark masses and the CKM elements are listed in Table 7. Values for the other observables are listed in Table 8. We do not assume unitarity of the CKM matrix for our analysis and our parameters are fit directly to the values determined by experimental measurements.
The relevant effective Hamiltonian for b → s + − is given by [50,51], where Here, = e, µ, τ . The Wilson coefficients induced by Z exchange are given by where i = 1, 2, 3 for = e, µ, τ , respectively. In Ref. [52] 4 , one or two of the Wilson coefficients are fitted to the latest data for R K ( * ) and b → s + − decay observables, while all the other Wilson coefficients are assumed to vanish. There are three scenarios in the one-dimensional analysis that have pulls larger than 5σ with respect to the SM prediction: For the two-dimensional analysis there are two patterns that have pulls larger than 6σ with respect to the SM prediction: In our analysis, we attempt to fit C ( ) =e,µ 9,10 to one of these 5 patterns. The central values and uncertainties of the other coefficients are assumed to be 0.0 ± 0.1. We also include imaginary parts, and try to fit them to 0.0±0.1. We remark here that non-zero imaginary parts are actually favored by the analysis of Ref. [61], however, we do not consider this possibility in the present paper.
C µ 9 is sizable in most of the above preferred patterns of Wilson coefficients. The Z contribution to C µ 9 can be estimated as where we note that (s 2 µ L + s 2 µ R ) ∼ 1 is required by a successful explanation of ∆a µ . Therefore, the b → s + − anomalies can be explained with small, O (10 −3 ), couplings of the SM quarks to the Z boson.
The Wilson coefficients with = τ contribute to the semi-leptonic decays B s → K ( * ) τ τ and B s → φτ τ . Branching fractions for these decays as a function of the Wilson coefficients are calculated in Ref. [69].

Neutral Meson Mixing
The neutral bosons, Z , χ and σ, give contributions to neutral meson mixing. We neglect contributions from the Z boson and the Higgs boson, since the flavor violating couplings are expected to be small, as discussed above and shown in Appendix B. We study K-K , B d -B d , B s -B s and D-D mixing [70][71][72][73][74][75].
The relevant effective Hamiltonian is given by where α, β are the color indices and The Wilson coefficients including O(α s ) corrections are given by [76] C VLL where µ is the MS-scheme renormalization scale. The gauge couplings g Z L,R and Yukawa couplings y S L,R are given by where the flavor indices are (i, j) = (2, 1), (3, 1), (3, 2) for K-K , B d -B d , or B s -B s mixing, respectively. For D-D mixing, (3.82) The right-right Wilson coefficients (C VRR ) are obtained by formally replacing L → R in the above expressions. The off-diagonal mixing matrix element in the respective meson's mass matrix is given by Here, m M is the mass of the meson M . Values of the operators at µ B = 1 TeV according to our own evaluation are listed in Table 9. For this we have used input values of meson masses, decay constants and quark masses which are listed in Table 10. Values of hadronic matrix elements are taken from the results of the respective lattice collaborations. We refer to Ref. [77] for hadronic matrix elements of K-K and f 2 BqB Bq . Those for f 2 Bq B (2-5) Bq and B (1-5) D are taken from Refs. [78,79] and Ref. [80], respectively. The QCD running between the respective lattice scales and µ = 1 TeV has been calculated based on the anomalous dimensions shown in Ref. [81]. SM contributions to K-K, B q -B q (q = d, s) mixing are given by Here, V is the 3 × 3 CKM matrix of the SM families. The Inami-Lim functions are given by (3.88) Short distance corrections are quantified by η 1,2,3 and η B . Values for all relevant factors and their respective references are listed in Table 10. The relevant observables are defined as Values of κ and φ are stated in Table 10. (∆M K ) exp is the experimentally determined value of the Kaon mass splitting and Γ D is the experimentally determined decay width of the D 0 meson; both are taken from the PDG [25].  [79,89] 163.53 ± 0.83 GeV The experiments measure the mass differences and K precisely. On the other hand, there are large theoretical uncertainties to estimate the SM contributions for these observables originating from the determination of the bag parameter, QCD factors and the CKM matrix elements. For K-K mixing, uncertainties come from η 1 , f K ,B K , κ and the CKM elements. The uncertainty of ∆M K is dominated by the NLO factor η 1 , while for K it is dominated by the CKM elements. With the Wolfenstein parametrization, K is approximately proportional to A 4 . Hence, we include the uncertainty from A = 0.836 ± 0.015 [25] as the CKM uncertainty together with those from f 2 KB K and κ . The relative uncertainties are estimated as 41% and 9.3% for ∆M K and K , respectively.
For the mass differences ∆M Bq ≡ ∆M q , we include the uncertainties originating from η B , f 2 BqB Bq and the absolute values of the CKM matrix elements. Note that unlike the analyses in Refs [75,90], we cannot reduce the uncertainties by assuming exact unitarity of the CKM matrix, because the unitarity of CKM matrix is not guaranteed in our model. Altogether, the relative uncertainties are estimated as 16% and 14% for ∆M d and ∆M s , respectively. For the CP asymmetry parameters S ψK S and S ψφ , we require our model to fit them within their experimental uncertainties.
For D-D mixing there is a large theoretical uncertainty from long-distance effects. The observed value is x D = 0.32 (14)% [44]. We simply require that the new physics contribution to x D should be less or equal than the size of the observed value, that is It is convenient to express the size of new physics contributions relative to the SM, , , (3.94) and these are given by Here,  Table 9 and neglecting the O (α s ) corrections in Eqs. (3.76)-(3.80). The SM contributions are calculated with the unitary CKM matrix fitted to the experimental values [25]. The coefficients for the lighter mesons tend to be larger because the SM contributions are smaller. Left-right contributions are enhanced, especially r K V LR , by the large hadronic matrix element itself and the enhancement by the running effects [81], see Table 9.
Similarly, ∆ NP x D is given by (3.103) We now comment on the box-diagram contributions involving W bosons and uptype quarks which are the dominant contribution in the SM. In general, the unitarity of the CKM matrix is violated by the mixing with the SU(2) L singlet VL quark. The GIM mechanism may, in principle, become invalid in our model. The mass independent contributions is proportional to a sum over the five internal quarks, This has the same structure as the weak-isospin part of the Z boson couplings. Using the analytical expressions of Appendix B, the size of flavor violating contribution is estimated as where D i is a mixing angle between the singlet VL quark D R and the SM down quark d i , and M D VL is a typical VL down quark mass. In addition, there can be an, in principle, important contribution which is enhanced by the heavy VL quark mass. Using Eq. (B.65), the dominant contribution is estimated as where U i is defined in the same way as D i . For B s -B s mixing, (i, j) = (3, 2), this should be compared with the top loop contribution λ Bs t 2 S 0 (x t ) ∼ 0.004 and is, thus, much smaller than the SM contribution. The same suppression happens for the other meson mixing. Thus, the violation of the GIM mechanism is extremely small.
The new bosons, in general, induce new physics contributions to B q → µ + µ − (q = d, s).
We refer to Ref. [48]. The relevant effective interactions are In this model, the coefficients are given by where i = 1, 2 for q = d, s, respectively. The SM contribution in C V AA is given by where η Y = 1.012 quantifies QCD corrections [91,92] and the loop function is given by (3.112) The decay width of B q → µ + µ − is proportional to We define the ratios of branching fractions of our model to the SM, Mind the bars: In the B s -B s system, the measured width difference between light and heavy mass eigenstates, y s := ∆Γ Bs / (2Γ Bs ) = 0.065 ± 0.005 [44], is not negligible [93,94]. The experimentally determined value for the branching ratio, therefore, corresponds to the time-integrated value where the mass-eigenstate rate asymmetry A ∆Γ is given by [95,96] Here, P s = |P s | e iφ P , S s = |S s | e iφ S and φ NP s relates to S ψφ , defined in Eq. (3.89), as 118) in the standard phase convention of the CKM matrix. A ∆Γ = 1 in the SM.
The SM predictions are [47,48], (3.120) The experimental values are [25], Altogether, the values of the ratios are given by The current data for BR (B s → µ + µ − ) has a slight tension with the SM prediction. We note that BR(B s → µ + µ − ) is included in the analysis of Ref. [52], where due to the tension a larger C µ 10 is favored. Nonetheless, we additionally include BR(B s → µ + µ − ) individually in our χ 2 analysis in order to take into account scalar contributions which were not included in [52].

Top Quark Decays
The mixing with the VL quarks may affect top quark decays. We study the dominant top quark decay t → W + b and the flavor violating decays t → Zq and t → hq (q = u, c).
The partial decay width and the branching fractions for the flavor violating decays are,

Z Physics
We now study potential signals of the Z gauge boson at the LHC, in gauge kinetic mixing, and in neutrino trident production. In general, there are exclusion regions from all these observables. Note that we do not include these observables in our χ 2 analysis, but only check a posteriori whether the respective constraints are fulfilled.

Dimuon Signals at the LHC
In the present model, the Z gauge boson should be lighter than about 800 GeV to explain ∆a µ . The most relevant Z -related process at the LHC is resonant dimuon production, (3.133) ∆a µ requires sizable couplings to muons, while small couplings to the SM quarks are enough to explain b → s + − anomalies. Hence, the Z boson will dominantly decay to muons and muon neutrinos, and its production cross section will be suppressed by the small couplings to the SM quarks. General LHC limits on Z bosons responsible for b → s + − anomalies are studied in Refs. [100,101]. Exclusion bounds are given in Ref. [

Gauge Kinetic Mixing
We assume that the gauge kinetic mixing between the U(1) Y and U(1) gauge boson is absent at tree-level. At the one-loop level mixing is unavoidable and the corresponding Z-Z mixing parameter is estimated as Current experimental limits are summarized in Ref. [106]. Values of ∼ 0.05 are not excluded provided that the Z is heavier than a few 100 GeV.

Neutrino Trident Production
The Z contributes to muon-neutrino induced muon pair production off a nucleus ν µ N → ν µ µ + µ − N , the so-called neutrino trident process [6,[107][108][109][110][111]. The cross section for this process at the CCFR experiment is estimated as [111] (see also [112] for a complete SM computation) The experimentally observed rate is σ CCFR /σ SM CCFR = 0.82 ± 0.28 at 95% C.L. The relevant effective interactions are where the neutrinos are taken as flavor states. In our model, the coupling constants are given by with Z boson contributions given by (3.138) Here, ĝ Z e L νµνµ is given by and we have used s 2 W = 0.23129 specifically for this process as in Ref. [111]. This constraint is relevant only for light Z 's and becomes unimportant for m Z 200 GeV.

χ 2 Fitting
We search for parameters that can explain both ∆a µ and b → s + − anomalies consistently with the other observables. For this, we attempt to minimize the χ 2 function, which are the Z mass, the VEV of φ, the U(1) gauge coupling constant, and the effective quartic couplings of the scalars Φ and φ. All other parameters are Yukawa coupling constants appearing in Eqs. (2.2)-(2.5). Generally, we assume that the Yukawa coupling constants are real, except for the couplings y u,d 13 , y u,d 31 which are taken to be complex in order to explain the complex phases in the CKM matrix. The Yukawa couplings involving the right-handed neutrinos with heavy Majorana masses, that is λ N i and y n ij , are not included in our χ 2 analysis as none of our 98 observables is sensitive to them. As discussed in Section 2.2.4, g < 0.35 is imposed, such that the gauge coupling stays perturbative up to ∼ 10 16 GeV. We restrict all Yukawa and effective quartic coupling constants to be smaller than unity and impose v φ ≤ 5.0 TeV.

Best Fit Points
We find a landscape of good fit points in similar phenomenological regions. We will focus our discussion on the four best fit points A, B, C and D with  are selected from points with the charged VL lepton heavier than 250 GeV and the fiducial cross section σ fid (pp → Z → µ + µ − ) smaller than the latest experimental limit. Point A is the global best fit point under these conditions. The point B is the best fit point of points with m E 1 > 1.2 TeV. This point has slightly larger χ 2 value than the other three best fit points (see Table 12), mainly because R th Bs→µµ ∼ 0.9 due to the smaller ReC µ

Phenomenology
We now discuss some global features of our model. Figure 3 shows fit points with χ 2 < 33 = N d.o.f. . The red circles (green triangles) have χ 2 < 28 (33). Points which are excluded by Z physics, namely LHC searches and/or neutrino trident production, are denoted by red crosses (green pluses) with the same color coding as above. The Z-Z mixing parameter is always less than or equal O (10 −3 ) and, thus, much smaller than the experimental limit. All subsequent plots show the same model parameter points as Figure 3.
The blue solid (dashed) line in Figure 3 corresponds to v Φ = 1.7 (2.0) TeV. Consistent with our analytical analysis of ∆a µ in Section 3.1.1, c.f. especially Eq. (3.18), there is no point with χ 2 < 28 (33) whenever v Φ > 1.7 (2.0) TeV. This results in an upper bound on the Z mass: m Z 840 GeV for g < 0. 35. We note that in Fig. 3 allowed and excluded points co-exist for similar values of m Z and g . This is because in this plane one does which have dramatic consequences for Z direct production as we will discuss now. Figure 4 shows the good fit points in the (m Z , σ fid (pp → Z → µµ)) plane, where we use the definition and cuts for the fiducial cross section σ fid of Ref. [102]. The blue solid line is the 95% C.L. limit from the ATLAS analysis [102]. Since the limit is given only for m Z > 250 GeV, we use an extrapolation down to lower masses shown by the dashed blue line. As a rough estimate for the sensitivity to be expected at the HL-LHC we can scale the limit on the cross section by 139/3000, the square root of the expected ratio of integrated luminosities. This sensitivity is shown as a purple, dot-dashed line in Fig. 4.

Z Physics
A small flavor violating coupling to Z , ĝ Z d L 23 ∼ 10 −3 is enough to explain the b → s + − anomalies. A diagonal coupling of Z to bottom quarks or to the light quarks could be sizable without changing other flavor violating observables. However, fitting the observed CKM matrix sets limits on the size of such couplings. Therefore, a good fit prefers small diagonal couplings to quarks. In agreement with that, our best fit points predict fiducial cross section roughly about an order of magnitude smaller than the current limits. We stress that the LHC limits were not part of the fit and only checked subsequently on good fit points.
Since ∆a µ requires sizable Z coupling to muons, a sizable muon neutrino coupling  is also predicted. Our model, therefore, is sensitive to the neutrino trident process if v Φ 350 GeV. Focusing on this mass range in Fig. 4, we see that there are a handful of points which are excluded exclusively by the trident constraints and not by LHC searches. On a different note, the one-loop induced gauge kinetic mixing for all points is O (10 −3 ) or less, much smaller than the current limits.

b → s + −
All the best fit points A-D are fitted to pattern (IV) ("C 9 and C 10 ", cf. eq. (3.68)). There are also a lot of points which are fitted to pattern (I) ("C 9 only", eq. (3.65)). We show our good fit points the (Re C µ 9 , Re C µ 10 ) plane in the left panel of Fig. 5. Points with pattern (IV) tend to have smaller χ 2 because of the tension in R Bs→µµ which favors non zero C 10 . The other patterns (II) ("C 10 only"), (III) ("C 9 = −C 10 "), and (V) ("C 9 and C 9 "), are hardly compatible with other observables and we will now discuss this in some detail. Making use of the analytic discussion in Appendix B, the Z couplings to muons can be expressed as  Figure 4: Good fit points in the (m Z , σ fid (pp → Z → µµ)) plane. The color coding is the same as in Fig. 3. The solid blue line is the 95% C.L. exclusion limit from ATLAS [102]. The blue dashed line extrapolates this line to test the points with m Z < 250 GeV. As purple dot-dashed line we show a rough estimate of the future sensitivity to be expected after HL-LHC. Excluded points below the blue line are not excluded by LHC direct searches but by neutrino trident production.
Hence, the ratio C µ 10 /C µ 9 is given by This indicates |C µ 10 /C µ 9 | ≤ 1, and that pattern (II) ("C 10 only") can never be realized. pattern (III) is C µ 9 ≈ −C µ 10 , implying s 2 µ R 1. However, as ∆a µ ∝ s µ L s µ R /v 2 Φ (cf. Eqs. (3.14), (3.15)) it would be suppressed in this case unless the suppression is compensated by a small v Φ 500 GeV. We show this on the right panel of Fig. 5, where one can clearly see that there are no good points with Re(C µ 10 /C µ 9 ) −0.8 for v Φ 1.0 TeV. Finally, pattern (V) is incompatible with neutral meson mixing: As can be seen from Eqs. (3.95), mixed LR contributions of Z exchange are enhanced by large negative coefficients. Since pattern (V) requires that ĝ Z d L 23 and ĝ Z d R 23 have opposite signs, their LR contribution to meson mixing adds constructively with the SM. As the SM prediction for ∆M s is already larger than the experimentally measured value, Z couplings compatible with pattern (V) would only ever increase the tension with experiment. This could be overcome if there were  Figure 5: Good fit points in the (−Re C µ 9 , Re C µ 10 ) and (v Φ , Re(C µ 10 /C µ 9 )) planes. The color coding is the same as in Fig. 3. The black line in the right plot correspond to pattern (IV) of Eq. (3.68). sizable negative contributions from the scalar exchange, however, the scalar couplings in our model are always suppressed as shown in Appendix B.

Standard Model Quark Sector
We fit our model parameters to match the quark masses, absolute values of the CKM matrix elements, and relative physical phases α, β and γ. We do not assume unitarity of the CKM matrix and fit our parameters directly to the experimentally determined absolute values and angles. In addition, we require our model to fit the Wilson coefficients of b → s + − processes such that the anomalies are matched. Furthermore, we fit to CPeven and CP-odd observables in K-K , B d -B d , B s -B s and D-D mixing as well as R νν K ( * ) , R B d(s) →µµ , BR(B s → Kτ τ ) and top quark decays. Of course, to some extent this approach consists of a "double fitting" as CKM angles and phases are themselves extracted also from some of these observables under the assumption of the SM. However, our approach should be valid here as NP contributions to the relevant observables are typically less than O(10%).
Our best fit values for CKM matrix elements and angles, relative to the SM extraction, are shown in Fig. 6. The brown lines and yellow bands show central values and their uncertainties as obtained in a global fit to the SM [25]. It is an important non-trivial crosscheck of our fitting procedure that we reproduce the SM best-fit values for most elements. In general, our results are consistent with the SM as most of the values agree within their 1σ region. However, some elements, namely |V cb |, |V td | and |V ts | show consistent deviations from the SM extraction. Perhaps these deviations could be tested by |V ud ||V us ||V ub ||V cd ||V cs ||V cb ||V td ||V ts ||V tb | α sin2β γ β s future experiments, which would be especially interesting if they are correlated with other observables. In Fig. 7 we show the good fit points in the (∆M d , ∆M s ) plane compared to experimental measurements and SM prediction with and without the assumption of CKM unitarity. Uncertainties of the SM predictions are shown in the figure. As discussed Section 3.3.2, the relative uncertainties for ∆M d and ∆M s are estimated as 16 % and 14 % without assuming unitarity of the CKM matrix. These uncertainties reduce to 9.8 % and 7.1 % if CKM unitarity is assumed [25].
Although there are sizable Z contributions these are still small compared with the uncertainties originating from the CKM elements without assuming unitarity. The CP asymmetry parameters S ψKs and S ψφ are well fit to the experimental values, cf. their values at the best fit points in Table 12 and experimental values in Table 8.
The right panel of Fig. 7 shows the good fit points in the (∆M K , K ) plane. Similar  to B-B mixing, the Z contributions are smaller than the uncertainties from the CKM values and QCD corrections.
There may be a sizable contribution from Z exchange to D-D mixing as well. However, theoretical errors here are too large to hope for a discrimination of SM and NP effects. Also the effects in top quark decays, including flavor violating ones, are too small to be tested by experiment.
Potentially, there are also contributions from the new scalar fields. However, as shown in Appendix B, the Yukawa coupling of the new scalars to the SM fermions first arise at the second order of the small mixing angles. Therefore, scalar contributions are very suppressed for B s → µµ. The ratio R Bs→µµ is predominantly determined by C µ 10 where a larger C µ 10 relaxes the tension between measurements and the SM prediction, see Table 12. Finally, BR(B s → Kτ τ ) is generally not much affected as the Z coupling to τ 's can be suppressed. At all the best fit points, R νν 1.06 which is a deviation much smaller than the experimental sensitivities.

Charged Lepton Flavor Violation
We now discuss predictions for charged LFV in this model. While we have included the experimental upper bounds on charged LFV in the fit, it still is interesting to see the size and spread of LFV obtained in this model. Figure 8 shows the best fit points in the (BR(µ → eγ), BR(µ → eee)) and (BR(τ → µγ), BR(τ → µµµ)) planes, respectively. For   Figure 8: Predictions of the best fit points for BR(µ → eγ), BR(µ → eee), BR(τ → µγ), as well as BR(τ → µµµ) and their correlations. The color coding is the same as in Fig. 3. The black lines are the current experimental 90% C.L. exclusion limits.
LFV muon decays, BR(µ → eγ) is much larger than BR(µ → eee). As already discussed in Section 3.1.3 and 3.1.4 this can be understood analytically. While the former decay is quadratically proportional to the tiny mixing angle e between electron and VL leptons, the latter decay scales with 6 e . Thus, our model predicts BR(µ → eγ) BR(µ → eee). For LFV τ decays, in contrast, BR(τ → µµµ) is roughly of the same order of magnitude as BR(τ → µγ). This can be understood because both of them are scaling as 2 τ , where τ is the small mixing angle between τ and the VL leptons. All the other LFV tau decays are suppressed by additional powers of τ and/or e . We see that especially for µ → eγ there are many best fit points close to the current upper limit. This limit will be significantly improved to < 6 × 10 −14 by the upcoming MEG II experiment [113]. Similarly an improvement on the upper limit of BR(τ → µγ) by up to two orders of magnitude is expected from Belle II [114]. Nonetheless, we find good fit points that extend into regions which will not be probed by upcoming experiments. We therefore conclude that while a future excess in the charged LFV channels µ → eγ and/or τ → µγ could consistently be explained in our model, those observables will not be the first to exclude this model.
Regarding charged LFV decays of SM bosons at tree level, we find that the respective branching fractions are more than several orders of magnitude smaller than the current bounds. In fact, the couplings of SM bosons to SM fermions are very close to their SM values which we have already discussed in Section 3.2 based on our analytical treatment shown in Appendix B.

Signals of Vector-Like Leptons
Finally, let us investigate collider signatures of VL fermions. As discussed in Section 3.1.1, the VL lepton mass is constrained to explain the muon anomalous magnetic moment, and ∆a µ = 2.68 × 10 −9 can be realized only within the parameter space shown in Fig. 2. In the same figure we show the masses of the lightest VL lepton m E 1 and m Z for our best fit points. Most points have m E 1 800 GeV and the points with m E 1 > 800 GeV are found only where m Z ∼ 500 GeV as expected from our analytical discussion illustrated by the contours in Fig. 2. In the upper panels of Fig. 9 we show the distribution of the heavier VL lepton masses m E 2 with respect to m Z (left) and with respect to the lighter VL lepton m E 1 (right). The black thick, thin, and dashed lines show mass splittings ∆m E := m E 2 −m E 1 = 0, 174, and 2·174 GeV, respectively. The mass splitting is bounded by ∼ λ e v H , and it is typically not very large, since the loop function C Z contributing to ∆a µ , defined in Eq. (3.16), is maximized when the masses are degenerate. Consequently, the heavier VL lepton E 2 is typically not much heavier than about 1.5 TeV. According to Ref. [115], the production cross section of a doublet VL lepton with mass 1.5 TeV is about O (0.1) fb which corresponds to about 30 (300) total events at the end of LHC (HL-LHC). Therefore, the HL-LHC could exclude the whole parameter space compatible with ∆a µ if the signals for VL leptons are very clean.
In the present model, high-multiplicity muon signals are expected from the production and decay of VL leptons. The decay modes crucially depend on the masses of the Z and χ bosons. We show the dominant two-body decay modes and their branching fractions at our best fit points in Tables 18-21 in Appendix D. If either of the following final states is kinematically allowed, the lightest VL lepton decays dominantly to E 1 → Z + µ, and/or E 1 → χ + µ. (4.5) For illustration, the lower left panel in Fig. 9 shows the sum of BR(E 1 → Z µ)+BR(E 1 → χµ) in dependence of m E 1 for our good fit points, cf. also Tables 18-21. Often either χ or Z is lighter than the VL leptons, as is the case for the best fit points A, B, and D. This comes about because a light scalar χ is favored to suppress its destructive contribution to ∆a µ , while the Z mass is controlled by the overall scale v Φ which is rarely above ∼ 1 TeV. On the contrary, if m E 1 < m Z , m χ (as in our best fit point C), E 1 decays predominantly to a SM boson and a muon or neutrino, The detailed ratio of these branching fractions depends on the mixing between the singletlike and doublet-like VL states. The final states in Eq. (4.6) have been studied as a signature of VL leptons in several papers [115][116][117][118][119][120][121]. However, there is no study by LHC experiments of VL leptons decaying to a muon based on LHC Run-2 data. A dedicated analysis of the experimental data shows that a singlet-like VL lepton at the point C with mass above 200 GeV can not be excluded [117].     Figure 9: Good fit points and their predicted values of the observables, m Z , m E 1 , m E 2 , m D 1 , and BR(E 1 → Z /χ + µ). The color coding is the same as in Fig. 3.
The final states in Eq. (4.5) have not been considered so far; they give rise to characteristic multi-lepton final states. This comes about because the Z boson predominantly decays to a pair of muons or muon neutrinos, cf. Tables 18-21. The scalar χ couples to SM fermions through the tiny Yukawa couplings induced by mixing effects, ∼ f SM m f SM /M F VL . However, at most points χ has a large coupling to muons because of the large mixing. In addition, couplings to third generation quarks could also be large due to an enhancement by their heavy masses. This is the case at our best fit point D, cf. Table 21. The sizable branching fractions of the exotic boson to pairs of muons provide clean resonance signals, Therefore, processes with dramatic multi-resonant multi-lepton final states, such as are expected from the VL lepton production. Here, (µ + µ − ) B is a pair of muons with a resonant feature at the invariant mass m 2 µµ ∼ m 2 B . At point D, the doublet-like VL neutrino also decays to the exotic scalar and signals such as are expected. This signal is expected if the lightest VL lepton is doublet-like. The heavier VL lepton decays in a more complicated way. It will predominantly decay to the lighter VL lepton under the emission of a SM boson, since there is sizable mixing between the VL leptons in order to enhance the left-right effects to ∆a µ . The most dramatic case may be which results in five SM leptons from one VL lepton. A pair of VL leptons, thus, could produce up to ten leptons per event. Although it may be difficult to reconstruct all of them, such high-multiplicity lepton signals provide a very distinctive event topology.

Signals of Vector-Like Quarks
The lower right panel of Fig. 9 shows the good fit points in the (m Z , m D 1 ) plane. Unlike for the VL leptons, there is no stringent upper limit on the VL quark masses. This is because small Z couplings to the SM quarks are enough to explain the b → s + − anomalies. Moreover, the mixing itself is given by and is independent of the Higgs VEV. The VL quarks can be within the reach of current and future collider experiments if they are light. For instance, point A has a singlet-like VL quark with mass ∼ 1.6 TeV. Since v Φ < 1.7 TeV is required, VL quark decays to Z or χ is always kinematically allowed. As for the VL leptons, dramatic signals involving Z or χ, e.g.
are expected. These high-multiplicity leptons with resonant features in association with jets provide another distinctive signal of this model.

Summary
We have presented an extension of the Standard Model with the addition of a vectorlike family of quarks and leptons which also carry a new U(1) charge. The model is constructed to address the known experimental anomalies associated with muons, i.e. the anomalous magnetic moment of the muon and the decays b → s + − . SM quarks and leptons feel the new Z gauge interactions only via mass mixing with the VL family. The model contains two additional SM singlet scalars, one that is used to model the masses of the VL family and another one that mixes the VL family with the SM states and obtains a VEV that spontaneously breaks the U(1) symmetry. We performed a global χ 2 analysis of the data, with 65 arbitrary model parameters and taking into account 98 observables. We have found many model points which satisfy χ 2 /N d.o.f. ≤ 1. We cannot simultaneously fit the anomalous magnetic moment of the electron and muon, because the theory is severely constrained by the experimental upper bound on BR(µ → eγ). 5 We, therefore, decided to only fit ∆a µ . All good fit points have BR(µ → eγ) > 10 −16 and BR(τ → µγ) > 10 −15 with the latter being strongly correlated with BR(τ → 3µ). Roughly half of our best fit points have BR(µ → eγ) in a range that is accessible by upcoming experiments. However, note that this is not a necessity of the model, i.e. BR(µ → eγ) could always be suppressed by tuning e L,R to zero without affecting the explanation of the anomalies or SM predictions. With regards to b → s + − decay processes, we fit the Wilson coefficients for new physics contributions as discussed in Ref. [52]. Only two of the five possible good fit points of this analysis can be fit in our model. The flavor violating decays of SM bosons, i.e. Higgs, W and Z are severely suppressed in our model. The vector-like quark induced coupling to Z also gives sizable contributions to neutral meson mixing, particularly K-K, B d,s -B d,s , and D-D mixing. The best-fit values for many CKM elements in our model consistently deviate from their experimental central values at the level of 1-2σ (as they do also in the Standard Model). Hence, more accurate constraints of CKM non-unitarity and more precise measurements of CKM elements would be very welcome to further test the model.
In order to understand the predictions for new physics we have presented four "best fit points" -A, B, C, and D with the masses of the new particles and their decay rates given in Tables 18, 19, 20 and 21, respectively. Many more details are given in the Appendices. The fit values for some selected observables are given in Table 12. Although the Z mass is typically significantly less than a TeV and it decays with a significant branching fraction to µ + µ − , we find many points not excluded by recent ATLAS searches for a dimuon resonance. We are also constrained by neutrino trident processes. The VL leptons are typically light, while the VL quarks are significantly heavier with mass of order a few TeV. Since the lightest VL leptons at best fit points, A, C and D, have mass between 300 -600 GeV, these states may be accessible even at the LHC, and even more so at the HL-LHC. They typically result in multi-muon production channels as discussed in Section 4.3.5.
This model is a prototype which highlights that fixing anomalies with consistent mod-els, while maintaining the successful Standard Model predictions, comes at a price: The model appears eminently testable and, therefore, can be excluded in many complementary ways.

Acknowledgments
We thank R.
Dermisek for useful discussions. The work of J.K. and S.R.

A Decay Width Formulas
Widths of two-body decays with both left-handed and right-handed interactions are summarized in this Appendix.

A.1 Scalar Decays
With the Yukawa interactions of a real scalar field φ and two fermions f 1 , f 2 , the partial width of φ is given by

A.2 Gauge Boson Decays
Gauge interactions of a vector field V , two fermions f 1 , f 2 , give the partial width,

A.3 Fermion Decays
If m f 2 > m f 1 + m φ , a fermion f 2 can decay as f 2 → f 1 φ through the Yukawa interaction in Eq. (A.1). The partial width is given by The partial width is given by where y = m 2 f 1 /m 2 f 2 and z = m 2 V /m 2 f 2 .

B Analytical Analysis
Many analytical formulae used in the main text are derived in this Appendix.

B.1 Diagonalization of the Dirac Mass Matrices
The 5 × 5 Dirac mass matrices are given by Here, the four-component vectors obey the following conditions, The vectors z L i , z R i can be any four-component vectors which satisfy Eq. (B.4). Another set of with arbitrary unitary matrices u L , u R also satisfy the conditions Eq. (B.4). We define these vectors such that the upper-left 3 × 3 block is diagonalized by using this degree of freedom. Rotating the mass matrix by these unitary matrices while keeping O (v H ) elements we obtaiñ where the 4 × 4 matrixm f is defined aŝ The matrixM f is 3 ⊕ 2 block diagonal, exceptỹ f L i ,R j v H . The mixing effects induced by these elements are O()ỹ f L,R v H /M F L,R ), suppressed by Yukawa couplings and VL fermion masses.
Next, we diagonalize the lower-right block. We are interested in parameters where M E L,R 250 GeV and M Q L,R 1.5 TeV to be consistent with LHC searches. The VL quarks are substantially heavier than v H , while the VL leptons are at a few hundred GeV with a Yukawa coupling λ e ∼ O (1) as required in order to explain ∆a µ . Fortunately, the other Yukawa couplings are small enough due to the small charged leptons masses.
Keeping λ e v H , the next order of unitary matrices is parametrized as with angles θ L,R that satisfy The rotated mass matrix is (B.10) We now give approximate forms for angles θ L,R and massesM F L,R . Here, we neglect the mixing angles and masses are approximately given by with an expansion parameter defined as Clearly, this approximation becomes inaccurate if the VL fermions are nearly massdegenerate. For the nearly mass-degenerate case one can proceed as follows. We introduce three mass parameters,

the masses are given bỹ
14) The mixing angles are given by Multiplying these unitary matrices one finds The second-order corrections to the upper-left block are given by Here, it does not matter whether Eq. (B.11) or Eq. (B.15) is used in the second equality; both give the same result at this accuracy. These corrections are negligibly small in the relevant parameter space. Altogether, the approximate mass basisf L ,f R is defined aŝ

B.2 EW Boson Couplings
Couplings of the fermions to SM bosons are completely SM-like at the leading order. The leading order unitary matrices of Eq. (B.2) transforms the SU(2) L gauge couplings as Here, V ij is a 3 × 3 unitary matrix, which is an identity matrix for the Z boson couplings where f = f . Since U 1 f L,R do not affect SM fermion couplings, only the mixing via U 2 induces flavor violating couplings of SM generations. Their size is estimated as and M F VL is the typical scale of VL fermions.
The Higgs boson couplings are aligned with the mass matrix by the rotation via U 0 f L,R , Hence, the Yukawa couplings to SM fermions are diagonal in the mass basis ifỹ f L,R are neglected. The mixing U 2 f L,R induces flavor violating couplings of size

B.3 Charged Leptons
For charged leptons, let us start in a basis in which the upper-left 3 × 3 block is diagonal, such that SM-LFV effects are induced only by λ L,E i . We can achieve this form by redefining e L i , e R i . Such a redefinition does not change the Z and Z couplings. The W boson couplings are changed, but this can be absorbed by a redefinition of the neutrinos. The Yukawa couplings to the scalars, namely Higgs boson and χ, are still aligned with the mass matrix. In our analysis, we assume that all parameters in the charged lepton sector are real.
There should be sizable mixing between the muon and VL leptons to explain ∆a µ , while mixing involving e or τ can be small to avoid LFV processes. We introduce mixing parameters involving muon, and those for electron and tau, We expect e , τ 1 in order to suppress the LFV processes. With these parameterizations, the leading order unitary matrices are given by The diagonal structure of the upper-left 3 × 3 block holds as far as e L,R , τ L,R 1. The large mixing with the muon and VL leptons indicate that λ e ∼ y e 2 ∼ m µ /v H to explain the muon mass without fine-tuning. The Yukawa couplings in the off-diagonal block are given byỹ Hence, the perturbative corrections to the off-diagonal elements are at most, Consequently, the basis defined in Eq. (B.22) is very close to the mass basis and flavor violating couplings of the charged leptons to the SM bosons are strongly suppressed.
Using the above results, the Z couplings to charged leptons are given by The mixing effects from U 2 e L,R are neglected. In the same approximation, the off-diagonal Yukawa couplings of χ to VL and SM fermions are given by However the mass of σ is not bounded by v Φ implying that its effects can be decoupled for large v φ . Using above results, the leading contribution to ∆a µ from Z and χ boson can be estimated by Eq. (3.14). Here we want to give more details on the combination of loop functions C LR appearing there. The relevant combination of loop functions is defined as where x L,R :=M 2 E L,R /m 2 Z , y L,R :=M 2 E L,R /m 2 χ . This is straightforwardly obtained from the sum of Eqs. (3.3) and (3.11), while using our approximations above. For large enough mass splitting, λ e v H |∆ E | we can use Eq. (B.11) to simplify this expression to On the other hand, if the VL mass splitting is small, |∆ E | M E , we obtain the following formula by using Eq. (B.15), Relatively light VL leptons with large Higgs Yukawa couplings (necessary to explain ∆a µ by chiral enhancement) may contribute significantly to h → γγ. The Higgs couplings to VL leptons are approximately given by where λ e is neglected. The amplitude of h → γγ from the VL lepton loop is given by

B.4 Quarks
Before starting the analytical discussion of the quark couplings, let us obtain an estimate for the necessary size of such couplings. Based on above results one finds that is required in order to explain ∆a µ . The Z contribution to C µ 9 can then be estimated as Thus, the b → s + − anomalies can be explained even with the permille Z couplings to quarks. Let us start the discussion of quark mass diagonalization in a basis where the upper-left blocks of the quark mass matrices have already been diagonalized, The quark couplings to the Higgs, Z, and Z bosons is the same as in the gauge basis, while the W boson couplings (2.29) are changed to Here, V ij ∼ V CKM is expected because of the small mixings between VL and SM quarks. In general, the mass matrices can be diagonalized exactly in the same way as for the charged leptons in the previous section, but all mixing angles can be taken to be small. We define the small mixing angles, With this parametrization, the unitary matrices are Here, we keep terms of O 2 and linear for the other angles. The Yukawa couplings in the off-diagonal elements of (B.6) are given bỹ The perturbative correction to the SM up quark mass matrix is Hence, the basis defined as in Eq. (B.22) is to a very good approximation the mass basis. Perturbative corrections for the down quarks are even smaller due to their lighter masses.
The Z boson couplings are estimated as Thus, the b → s + − anomaly requires Q 2 Q 3 ∼ 10 −3 .
We now estimate the couplings to Z and Higgs boson. We focus on the up quark sector, where mixing effects are larger due to the heavy top quark. The Higgs Yukawa coupling are estimated as and the weak-isospin part of the Z boson couplings are where M U VL is a typical VL up quark mass. Thus, the EW boson couplings to the SM up quarks do not significantly deviated from their SM values. Those for the down quarks are even more suppressed by the smaller down quark masses.
The W boson couplings to the right-handed SM quarks are estimated as An estimation of the non-unitarity of the extended CKM matrix has already been given in Eq. (3.105). In addition, the off-diagonal elements involving the VL quarks are where we have neglected terms in U 1 u L,R , U 1 d L,R of O s u L,R , s d L,R . These effects are much smaller for quarks as compared to charged leptons owing to the heavier VL quark masses.

B.5 Neutrinos
We consider the type-I seesaw mechanism with Majorana masses M Maj ∼ 10 14 GeV.
with a unitary matrix V n L which is given by Here, whereM L is defined as the normalization factor of n † N n N = 1. The Dirac matrix after this rotation is given bỹ where a, b = 4, 5. Owing to the first rotation U 0 n L , the first three columns in the 4th and 5th rows inM n are all zero and there is no mixing between the singlet-like VL neutrinos (N R , N R ) and left-handed light neutrinosν L .
Back in the full 10 × 10 matrix (recall Eqs. (2.15) and (2.16))), we note that only the first three rows couple to ν R i and they will have O (M Maj ) masses. Therefore,ν L i are approximately the light Majorana neutrinos, whileÑ L ,Ñ L together with N R , N R form approximately Dirac neutrinos. After integrating out the O (M Maj ) right-handed neutrinos, the effective 7 × 7 mass matrix is given by where 2×3 , 2×2 , 3×2 are mass parameters suppressed by the Majorana mass term such as (µ n ) T M −1 Maj µ n . These terms are at most of O(v 2 Φ /M Maj ) and, therefore, negligible compared to M D . The 3 × 3 effective Majorana mass matrix m ν for the active neutrinosν L is given by Altogether, the mass terms for the neutrinos are decomposed as where the family indices are omitted. The mixing terms among the three types of neutrinos, namelyν L , ν R and the Dirac neutrinos, are all suppressed by M Maj and not stated. The remaining mass matrices then are diagonalized as The mass basis for the active neutrinosν L , heavy right-handed neutrinosν R and Dirac neutrinosN a := N a L ,N a R are given by where i, j = 1, 2, 3 , a, b = 4, 5 , and A = 1, 2, 3, 4, 5. Altogether, the relation between the mass and gauge basis is The mixing between the left-and right-handed neutrinos is suppressed by the heavy Majorana masses, and therefore negligible. The W boson coupling to the SM leptons is given by where k = 1, 2, 3. Therefore, we define flavor neutrino states ν f and the 3 × 3 PMNS matrix as . This would be an interesting possibility but is beyond the scope of the present paper. We restrict ourselves to the case where λ n is as small as λ e and the mixing between active and new Dirac neutrinos is negligible. In this case, the Z boson couplings to the flavor neutrinos defined in Eq. (B.81) is given by Hence, the coupling of Z to muon neutrinos can be estimated as where we have omitted couplings involving the Dirac neutrinos. The electron and tau neutrinos have tiny couplings to Z due to the tiny mixing angles e L , τ L . We can find parameters consistent with the experimental results on neutrino mass differences and PMNS mixing [25]. This can be done by fitting the Yukawa couplings involving the three generations of active neutrinos and the Majorana masses. The explicit values of the neutrino Yukawa couplings at the best fit points are shown in the subsequent section of the appendix. For simplicity, we always assume λ N i = 0 i , M ij Maj = M 0 δ ij and take the lightest neutrino to be massless. However, nothing in our analysis really depends on these assumptions. All of our fit points realize the experimental values of the mass differences,

D Full List of Observables at Best Fit Points
The CKM matrices at the best fit points are as follows: