Global Fits of Third Family Hypercharge Models to Neutral Current B-Anomalies and Electroweak Precision Observables

While it is known that third family hypercharge models can explain the neutral current $B-$anomalies, it was hitherto unclear whether the $Z-Z^\prime$ mixing predicted by such models could simultaneously fit electroweak precision observables. Here, we perform global fits of several third family hypercharge models to a combination of electroweak data and those data pertinent to the neutral current $B-$anomalies. While the Standard Model is in tension with this combined data set with a $p-$value of $.00068$, simple versions of the models (fitting two additional parameters each) provide much improved fits. The original Third Family Hypercharge Model, for example, has a $p-$value of $.065$, i.e. $\sqrt{\Delta \chi^2}=6.5\sigma$.

The current paper is about the third family hypercharge option. The Third Family Hypercharge Model [41] (henceforth abbreviated as the 'Y 3 model') explains the hierarchical heaviness of the third family and the smallness of quark mixing. It was shown to successfully fit NCBAs, along with constraints from B s −B s mixing and LFU constraints on Z 0 boson interactions. The ATLAS experiment at the LHC has searched 1 for the production of BSM resonances that yield a peak in the di-muon invariant mass (m µµ ) spectrum, but have yet to find a significant one [64]. This implies a lower bound upon the mass of the Z in the Y 3 model, M X > 1.2 TeV [65], but plenty of viable parameter space remains which successfully explains the NCBAs. A variant, the Deformed Third Family Hypercharge Model (DY 3 model), was subsequently introduced [43] in order to remedy a somewhat ugly feature (ugly from a naturalness point of view) in the construction of the original Y 3 model: namely, that a Yukawa coupling allowed at the renormalisable level was assumed to be tiny in order to agree with strict lepton flavour violation constraints. The DY 3 model can simultaneously fit the NCBAs and be consistent with the ATLAS di-muon direct search constraint for 1.2 TeV< M X <12 TeV. We will also present results for a third variant, the DY 3 model, which is identical to the DY 3 model but with charges for second and third family leptons interchanged. As we show later on, this results in a better fit to data due to the different helicity structure of the couplings of the Z boson to muons (see section 4.3 for details).
In either of these third family hypercharge models, the local gauge symmetry of the SM is 2 extended to SU (3) × SU (2) × U (1) Y × U (1) X . This is spontaneously broken to the SM gauge group by the non-zero vacuum expectation value (VEV) of a SM-singlet 'flavon' field θ that has a non-zero U (1) X charge. In each model, the third family quarks' U (1) X charges are equal to their hypercharges whereas the first two family quarks are chargeless under U (1) X . We must (since it is experimentally determined to be O(1) and is therefore inconsistent with a suppressed, non-renormalisable coupling) ensure that a renormalisable top Yukawa coupling is allowed by U (1) X ; this implies that the SM Higgs doublet field should have U (1) X charge equal to its hypercharge. Consequently, when the Higgs doublet acquires a VEV to break the electroweak symmetry, this gives rise to Z 0 − Z mixing [41]. Such mixing is subject to stringent constraints from electroweak precision observables (EWPOs), in particular from the ρ-parameter, which encodes the ratio of the masses of the Z 0 boson and the W boson [66].
Third family hypercharge models can fit the NCBAs for a range of the ratio of the Z gauge coupling to its mass g X /M X which does not contain zero. This means that it is not possible to 'tune the Z − Z mixing away' if one wishes the model to fit the NCBAs. As a consequence, it is not clear whether the EWPOs will strongly preclude the (D)Y 3 models from explaining the NCBAs or not.
The purpose of this paper is to perform a global fit to a combined set of electroweak and NCBA-type data, along with other relevant constraints on flavour changing neutral currents (FCNCs). It is clear that the SM provides a poor fit to this combined set, as Table 1 shows. A p−value 3 of .00068 corresponds to 'tension at the 3.4σ level'. The (D)Y 3 models are of particular interest as plausible models of new physics if they fit the data significantly better than the SM, a question which can best be settled by performing appropriate global fits.
Our paper proceeds as follows: we introduce the models and define their parameter spaces in §2. At renormalisation scales at or below M X but above M W , we encode the new physics effects in each model via the Standard Model Effective Theory (SMEFT). We calculate the leading (dimension-6) SMEFT operators predicted by our models at the scale M X in §3. These provide the input to the calculation of observables by smelli-2.2.0 [67] 4 , which we describe at the beginning of §4. The results of 2 Possible quotients of the gauge group by discrete subgroups play no role in our argument and so we shall henceforward ignore them. 3 All χ 2 and p−values which we present here are calculated in smelli-2.2.0. We estimate that the numerical uncertainty in the smelli-2.2.0 calculation of a global χ 2 value is ±1 and the resulting uncertainty in the second significant figure of any global p−value quoted is ±3. 4 We use the development version of smelli-2.2.0 that was released on github on 8 th March 2021 which we have updated to take into account 2021 LHCb measurements of RK , BR(Bs → µ + µ − ) and BR(B → µ + µ − ).  Table 1. SM goodness of fit for the different data sets we consider, as calculated by smelli-2.2.0. We display the total Pearson's chi-squared χ 2 for each data set along with the number of observables n and the data set's p−value. The set named 'quarks' contains BR(Bs → φµ + µ − ), BR(Bs → µ + µ − ), ∆ms and various differential distributions in B → K ( * ) µ + µ − decays among others, whereas 'LFU FCNCs' contains R K ( * ) and some B meson decay branching ratios into di-taus. Our sets are identical to those defined by smelli-2.2.0 and we refer the curious reader to its manual [67], where the observables are enumerated. We have updated RK and BR(B s,d → µ + µ − ) with the latest LHCb measurements as detailed in the text.
the fits are presented in the remainder of §4, before a discussion in §5.

Models
In this section we review the models of interest to this study, in sufficient detail so as to proceed with the calculation of the SMEFT Wilson coefficients (WCs) in the following section. Under SU (3)×SU (2)×U (1) Y , we define the fermionic fields such that they transform in the following representations: 3} is a family index ordered by increasing mass. Implicit in the definition of these fields is that we have performed a flavour rotation so that d Li , e Li , e Ri , d Ri , u Ri are mass basis fields. In what follows, we denote 3-component column vectors in family space with bold font, for example u L := (u L , c L , t L ) T . The Higgs doublet is a complex scalar φ ∼ (1, 2, +1/2), and all three models which we consider (the Y 3 model, the DY 3 model and the DY 3 model) incorporate a complex scalar flavon with SM quantum numbers θ ∼ (1, 1, 0), which has a U (1) X charge X θ = 0 and is used to Higgs the U (1) X symmetry, such that its gauge boson acquires a mass at the TeV scale or higher. In the following we will present our results for three variants of third family hypercharge models, which differ in the charge assignment for the SM fields: -The Y 3 model, introduced in [41]. Only third generation fermions have non-zero U (1) X charges. The charge assignments can be read in Table 2.
-The DY 3 model as introduced in [43]. It differs from the Y 3 model in that charges have been assigned to the second generation leptons as well, while still being anomaly free. -The DY 3 model, which differs from the DY 3 model in that the charges for third and second generation lefthanded leptons are interchanged. The charge assignments can be read in Table 3.
All three of these gauge symmetries have identical couplings to quarks, coupling only to the third family via hypercharge quantum numbers. This choice means that, of the quark Yukawa couplings, only the top and bottom Yukawa couplings are present at the renormalisable level. Of course, the light quarks are not massless in reality; their masses, as well as the small quark mixing angles, must be encoded in higher-dimensional operators that come from a further layer of heavy physics, such as a suite of heavy vector-like fermions at a mass scale Λ > M X /g X , where g X is the U (1) X gauge coupling. Whatever this heavy physics might be, the structure of the light quark Yukawa couplings will be governed by the size of parameters that break the U (2) 3 global := U (2) q × U (2) u × U (2) d accidental global symmetry [68][69][70][71][72] of the renormalisable third family hypercharge lagrangians. For example, a minimal set of spurions charged under both U (2) 3 global and U (1) X was considered 5 in Ref. [73], which reproduces the observed hierarchies in quark masses and mixing angles when the scale Λ of new physics is a factor of 15 or so larger than M X /g X . Taking this hierarchy of scales as a general guide, and observing that the global fits to electroweak and flavour data that we perform in this paper prefer M X /g X ≈ 10 TeV, we expect the new physics scale to be around Λ ≈ 150 TeV. This scale is high enough to suppress most contributions of the heavy physics, about which we remain agnostic, to low energy phenomenology including precise flavour bounds 6 . For this reason, we feel safe in neglecting the contributions of the Λ scale physics to the SMEFT coefficients that we calculate in § 3, and shall not consider it in any further detail.
Continuing, we will first detail the scalar sector, which is common to (and identical in) all of the third family hypercharge models, before going on to discuss aspects of each model that are different (most importantly, the couplings to leptons).

The scalar sector
The coupling of the flavon to the U (1) X gauge field is encoded in the covariant derivative where X µ is the U (1) X gauge boson in the unbroken phase and g X is its gauge coupling. The flavon θ is assumed to acquire a VEV θ at (or above) the TeV scale, which spontaneously breaks U (1) X . Expanding θ = ( θ + ϑ)/ √ 2, its kinetic terms (D µ θ) † D µ θ in the Lagrangian density give the gauge boson a mass M X = X θ g X θ through the Higgs mechanism. After electroweak symmetry breaking, the electrically-neutral gauge bosons X, W 3 and B mix, giving rise to γ, Z 0 and Z as the physical mass eigenstates [41]. To terms of order (M 2 Z /M 2 X ), the mass and the couplings of the X boson are equivalent to those of the Z boson. Because we take M X M Z , the matching to the SMEFT ( §3) should be done in the unbroken electroweak phase, where it is the X boson that is properly integrated out. In the rest of this section we therefore specify the U (1) X sector via the X boson and its couplings. Throughout this paper, we entreat the reader to bear in mind that in terms of searches and several other aspects of their phenomenology, to a decent approximation the X boson and the Z boson are synonymous.
The covariant derivative of the Higgs doublet is where W a µ (a = 1, 2, 3) are unbroken SU (2) gauge bosons, σ a are the Pauli matrices, g is the SU (2) gauge coupling, B µ is the hypercharge gauge boson and g is the hypercharge gauge coupling. The kinetic term for the Higgs field, (D µ φ) † (D µ φ), contains terms both linear and quadratic in X µ . It is the linear terms that, upon integrating out the X µ boson, will give the leading contribution to the SMEFT in the form of dimension-6 operators involving the Higgs, as we describe in §3.
The charges of the fermion fields differ between the Y 3 model and the DY 3 model, as follows.

Fermion couplings: the Y 3 model
The Y 3 model has fermion charges as listed in Table 2 (in the gauge eigenbasis), leading to the following Lagrangian density describing the X boson-SM fermion couplings [41]: where Λ (I) are Hermitian 3-by-3 matrices.
and Ω and Ψ are described in §2.3. The V I are 3-by-3 unitary matrices describing the mixing between fermionic gauge eigenstates and their mass eigenstates. Note that the quark doublets have been family rotated so that the d Li (but not the u Li fields) correspond to their mass eigenstates. Similarly, we have rotated L i such that e Li align with the charged lepton mass eigenstates, but ν Li are not. This will simplify the matching to the SMEFT operators that we perform in §3. We now go on to cover the X boson couplings in the DY 3 model before detailing the fermion mixing ansatz (which is common to all three models).

Fermion couplings: the DY 3 model
For the DY 3 model with the charge assignments listed in Table 3, the Lagrangian contains the following X boson-SM fermion couplings [43]: The matrices Λ (I)

Fermion mixing ansatz
The CKM matrix and the PMNS matrix are predicted to be respectively. For all of the third family hypercharge models that we address here, the matrix element (V (d L ) ) 23 must be non-zero in order to obtain new physics contributions of the sort required to explain the NCBAs. Moreover, in the Y 3 model we need (V e L ) 23 = 0 in order to generate a coupling (here left-handed) to muons 7 . These will lead to a BSM contribution to a Lagrangian density term in the weak effective theory proportional to (bγ µ P L s)(μγ µ P L µ), where P L is the left-handed projection operator, which previous fits to the weak effective theory indicate is essential in order to fit the NCBAs [17][18][19][20][21][22][23].
In order to investigate the model further phenomenologically, we must assume a particular ansatz for the unitary fermion mixing matrices V I . Here, for V d L , we choose the 'standard parameterisation' often used for the CKM matrix [74]. This is a parameterisation of a family of unitary 3 by 3 matrices that depends only upon one complex phase and three mixing angles (a more general parameterisation would also depend upon five additional complex phases): (10) where s ij := sin θ ij and c ij := cos θ ij , for angles θ ij , δ ∈ R/(2πZ). To define our particular ansatz, we choose angles such that (V d L ) ij = V ij for all ij = 23, i.e. we insert the current world-averaged measured central values of θ ij and δ [74], except for the crucial mixing angle θ 23 , upon which the NCBAs sensitively depend. Thus, we fix the angles and phase such that s 12 = 0.22650, s 13 = 0.00361 and δ = 1.196 but allow θ 23 to float as a free parameter in our global fits. Following Refs. [41,43], we choose simple forms for the other mixing matrices which are likely to evade strict FCNC bounds. Specifically, we choose V d R = 1, V u R = 1 and in the Y 3 model 8 , and V e L = 1 in the DY 3 model. Finally, V u L and V ν L are then fixed by (9) and the measured CKM/PMNS matrix entries. For the remainder of this paper, when referring to the Y 3 model or the DY 3 model, we shall implicitly refer to the versions given by this mixing ansatz (which, we emphasise, is taken to be the same for all third family hypercharge models aside from the assignment of V e L ). Next, we turn to calculating the complete set of dimension-6 WCs in the SMEFT that result from integrating the X boson out of the theory.

SMEFT Coefficients
So far, particle physicists have found scant direct evidence of new physics below the TeV scale. This motivates the study of BSM models whose new degrees of freedom reside at the TeV scale or higher. In such scenarios, it makes sense to consider the Standard Model as an effective field theory realisation of the underlying high energy model. If one wishes to remain agnostic about further details of the high energy theory, this amounts to including all possible operators consistent with the SM gauge symmetries and performing an expansion in powers of the ratio of the electroweak and new physics scales.
The Standard Model Effective Field Theory [75][76][77] is such a parameterisation of the effects of heavy fields beyond the SM (such as a heavy X boson field of interest to us here) through d > 4 operators built out of the SM fields. In this paper we will work with operators up to dimension 6 (i.e. we will go to second order in the power expansion). Expanding the SMEFT to this order gives us a very good approximation to all of the observables we consider; the relevant expansion parameter for the EWPOs is (M Z /M X ) 2 1, and in the case of observables derived from the decay of a meson of mass m, the relevant expansion parameter is (m/M X ) 2 1. By restricting to the M X > 1 TeV region, we ensure that both of these mass ratios are small enough to yield a good approximation.
The SMEFT Lagrangian can be expanded as where O 5 schematically indicates Weinberg operators with various flavour indices, which result in neutrino masses and may be obtained by adding heavy gauge singlet chiral fermions to play the role of right-handed neutrinos. The sum that is explicitly notated then runs over all mass dimension-6 SMEFT operators, and the ellipsis represents terms which are of mass dimension (in the fields) greater than 6. The WCs C i have units of [mass] −2 . In the following we shall work in the Warsaw basis, which defines a basis in terms of a set of independent baryon-numberconserving operators [78]. By performing the matching between our models and the SMEFT, we shall obtain the set of WCs C i at the scale M X , which can then be used to calculate predictions for observables. To see where these dimension-6 operators come from, let us first consider the origin of four-fermion operators. We may write the fermionic couplings of the underlying theory of the X boson, given in (4) and (7) for the Y 3 model and DY 3 model respectively, as where is the fermionic current that the X boson couples to, where the sum runs over all pairs of SM Weyl fermions ψ i . The couplings κ ij are identified from (4) or (7), depending upon the model. After integrating out the X boson in processes such as the one in Fig. 1, one obtains the following terms in the effective Lagrangian: We match the terms thus obtained with the four-fermion operators in the Warsaw basis [78] in order to identify the 4-fermion SMEFT WCs in that basis. These 4-fermion operators are not the only SMEFT operators that are produced at dimension-6 by integrating out the X bosons of our models. Due to the tree-level U (1) X charge of the SM Higgs, there are also various operators in the Higgs sector of the SMEFT, as follows. The (linear) couplings of the X boson to the Higgs, as recorded in (3), can again be written as the coupling of X µ to a current, viz.
where this time is the bosonic current to which the X boson couples, where Due to the presence of X boson couplings to operators which are bi-linear in both the fermion fields (J µ ψ ) and the Higgs field (J µ φ ), integrating out the X bosons gives rise to cross-terms which encode dimension-6 operators involving two Higgs fields, one SM gauge boson, and a fermion bi-linear current. Diagrammatically, these operators are generated by integrating out the X boson from Feynman diagrams such as that depicted in Fig. 2.
Finally, there are terms that are quadratic in the bosonic current J µ φ , which encode dimension-6 operators involving four Higgs fields and two SM covariant derivatives. The corresponding Feynman diagram is given in Fig. 3. This accounts, schematically, for the complete set of dimension-6 WCs generated by either the Y 3 model or the DY 3 model 9 . We tabulate all the non-zero WCs generated in this way in Table 4 for the Y 3 model and in Table 5 for the DY 3 model.

Global Fits
Given the complete sets of dimension-6 WCs (Tables 4  and 5) as inputs at the renormalisation scale M X , 10 we use the smelli-2.2.0 program to calculate hundreds of observables and the resulting likelihoods. The smelli-2.2.0 program is based upon the observable calculator flavio-2.2.0 [79], using Wilson-2.1 [80] for running and matching WCs using the WCxf exchange format [81].
In a particular third family hypercharge model, for given values of our three input parameters θ 23 , g X , M X , the WCs in the tables are converted to the non-redundant 9 Note that in deriving the WCs we have assumed that the kinetic mixing between the X boson field strength and the hypercharge field strength is negligible at the scale of its derivation, i.e. at MX . 10 Strictly speaking, the parameters gX and θ23 that we quote are also implicitly evaluated at a renormalisation scale of MX throughout this paper.

WC
value WC value Table 4. Non-zero dimension-6 SMEFT WCs predicted by the Y3 model, in units of g 2 X /M 2 X , in the Warsaw basis [78]. We have highlighted the coefficient (for i = 2, j = 3) that is primarily responsible for the NCBAs in bold font.

WC
value WC value  Table 5. Non-zero dimension-6 SMEFT WCs predicted by the DY 3 model, in units of g 2 X /M 2 X , in the Warsaw basis [78]. We have highlighted the coefficient that is primarily responsible (for i = 2, j = 3) for the NCBAs in bold font. WCs for the original DY3 model may be obtained by switching the l indices 2 ↔ 3 everywhere. We have updated the data used by flavio-2.2.0 with 2021 LHCb measurements of BR(B d,s → µ + µ − ) taken on 9 fb −1 of LHC Run II data [8] by using the two dimensional Gaussian fit to current CMS, ATLAS and LHCb measurements presented in Ref. [83]: with an error correlation of ρ = −0.27. The most recent measurement by LHCb in the di-lepton invariant mass squared bin 1.1 < Q 2 /GeV 2 < 6.0 is where the first error is statistical and the second systematic [3] (this measurement alone has a 3.1σ tension with the SM prediction of 1.00). We incorporate this new measurement by fitting the log likelihood function presented in Ref. [3] with a quartic polynomial. The SMEFT weak scale WCs are then matched to the weak effective theory and renormalised down to the scale of bottom mesons using QCD×QED. Observables relevant to the NCBAs are calculated at this scale. smelli-2.2.0 then organises the calculation of the χ 2 statistic to quantify a distance (squared) between the theoretical prediction and experimental observables in units of the uncertainty. In calculating the χ 2 value, experimental correlations between different observables are parameterised and taken into account. Theoretical uncertainties are modelled as being multi-variate Gaussians; they include the effects of varying nuisance parameters and are approximated to be independent of new physics. Theory uncertainties and experimental uncertainties are then combined in quadrature.
We note that an important constraint on Z models that explain the NCBAs is that from ∆m s (included by 11 Jupyter notebooks (from which all data files and plots may be generated) have been stored in the anc/ subdirectory of the arXiv version of this paper.  Table 6. Goodness of fit for the different data sets we consider for the Y3 model as calculated by smelli-2.2.0 for MX = 3 TeV. We display the total χ 2 for each data set along with the number of observables n and the data set's p−value.
The data set named 'quarks' contains BR(Bs → φµ + µ − ), BR(Bs → µ + µ − ), ∆ms and various differential distributions in B → K ( * ) µ + µ − decays among others, whereas 'LFU FC-NCs' contains R K ( * ) and some B meson decay branching ratios into di-taus. Our data sets are identical to those defined by smelli-2.2.0 and we refer the curious reader to its manual [67], where the observables are enumerated. We have updated to RK and BR(B s,d → µ + µ − ) with the latest LHCb measurements as detailed in §4. Two free parameters of the model were fitted: θ23 = −0.145 and gX = 0.426.
smelli-2.2.0 in the category of 'quarks' observables), deriving from measurements of B s − B s mixing, because of the tree-level BSM contribution to the process depicted in Fig. 4. The impact of this constraint has significantly varied over the last decade, to a large degree because of numerically rather different lattice or theory inputs used to extract the measurement [84][85][86]. Here, we are wedded to the calculation and inputs used by smelli-2.2.0, allowing some tension in ∆m s to be traded against tension present in the NCBAs. As we shall see, in all the models that we consider the global fit is fairly insensitive to M X , provided we specify M X > 2 TeV or so in order to be sure to not contravene ATLAS di-muon searches [43,64,65]. We will demonstrate this insensitivity to M X below (see Figs. 9 and 15), but for now we shall pick M X = 3 TeV and scan over the pair (g X × 3 TeV/M X ) and θ 23 . Since the WCs at M X all scale like g X /M X , the results will approximately hold at different values of M X provided that g X is scaled linearly with M X . The running between M X and the weak scale breaks this scaling, but such effects derive from loop corrections ∝ (1/16π 2 ) ln(M X /M Z ) and are thus negligible to a good approximation.

Y 3 model fit results
The result of fitting θ 23 and g X for M X = 3 TeV is shown in Table 6 for the Y 3 model. The 'global' p−value is calculated by assuming a χ 2 distribution with n − 2 degrees of freedom, since two parameters were optimised. The fit is encouragingly of a much better quality than the one of the SM. We see that the fits to the EWPOs and NCBAs are simultaneously reasonable.
The EWPOs are shown in more detail in Fig. 5, in which we compare some pulls in the SM fit versus the Y 3 model best-fit point. We see that there is some improvement in the prediction of the W -boson mass, which the SM fit predicts is almost 2σ too low (as manifest in the ρ-parameter being measured to be slightly larger than one [74], for M Z taken to be fixed to its SM value). The easing of this tension in M W is due precisely to the Z − Z mixing in the (D)Y 3 models. The non-zero value of the SMEFT coefficient C φD breaks custodial symmetry, resulting in a shift of the ρ-parameter away from its treelevel SM value of one, to [66] where v is the SM Higgs VEV. Rather than being dangerous, as might reasonably have been guessed, it turns out that this BSM contribution to ρ 0 is in large part responsible for the Y 3 model fitting the EWPOs approximately as well as the SM does.
We also see that σ had 0 , the e + e − scattering crosssection to hadrons at a centre-of-mass energy of M Z , is better fit by the Y 3 model than the SM. Although the other EWPOs have some small deviations from their SM fits, the overall picture is that the Y 3 model best-fit point has an electroweak fit similar to that of the SM. Shaded regions are those preferred by the data set in the legend at the 95% confidence level (CL). The global fit is shown by the solid curves, where the inner (outer) curves show the 70%(95%) CL regions, respectively. The set named 'quarks' contains BR(Bs → φµ + µ − ), BR(Bs → µ + µ − ), ∆ms and various differential distributions in B → K ( * ) µ + µ − decays among others, whereas 'LFU FCNCs' contains R K ( * ) and some B meson decay branching ratios into di-taus. Our sets are identical to those defined by smelli-2.2.0 and we refer the curious reader to its manual [67], where the observables are enumerated. We have updated RK and BR(B s,d → µ + µ − ) with the latest LHCb measurements as detailed in §4. The black dot marks the locus of the best-fit point. In order to see which areas of parameter space are favoured by the different sets of constraints, we provide Fig. 6. The figure shows that the EWPOs and different sets of NCBA data all overlap at the 95% CL. The bestfit point has a total χ 2 of 43 less than that of the SM and is marked by a black dot. The separate data set contributions to χ 2 at this point are listed in Table 6. In order to calculate 70% (95%) CL bounds in the 2-dimensional parameter plane, we draw contours of χ 2 equal to the best-fit value plus 2.41 (5.99) respectively, using the combined χ 2 incorporating all the datasets.
We further study the different χ 2 contributions for the   9. From Fig. 7, we see that large couplings g X > 0.6 are disfavoured by EWPOs as well as the NCBAs. From Fig. 8 we see that the EWPOs are insensitive to the value of θ 23 in the vicinity of the best-fit point but the NCBAs are not. At large −θ 23 the Y 3 model suffers due to a bad fit to the B s −B s mixing observable ∆m s . In Fig. 9, we demonstrate the approximate insensitivity of χ 2 near the best-fit point to M X , provided that g X is scaled linearly with it.
Finally, we display some individual observables of interest in Fig. 10 at the Y 3 model best-fit point. While some of the prominent NCBA measurements (for example R K in the bin of m 2 µµ between 1.1 GeV 2 and 6 GeV 2 ) fit considerably better than the SM, we see that this is partly compensated by a worse fit in ∆m s , as is the case for many Z models for the NCBAs. The P 5 observable (derived in terms of angular distributions of B 0 → K * µ + µ − decays [87,88]) shows no significant change from the SM prediction in the bin that deviates the most significantly from experiment: m 2 µµ ∈ (4, 6) GeV 2 , as measured by LHCb [12] and ATLAS [13]. The fit to BR(Λ b → Λµ + µ − ) is slightly worse than that of the SM in one particular bin, as shown in the figure. Some other flavour observables in the flavour sector, notably various bins of BR(B → K ( * ) µ + µ − ), show some small differences in pulls between the SM and the Y 3 model. Whilst there are many of these and in aggregate they make a difference to the overall χ 2 , there is no small set of observables that pro-   Table 7. Goodness of fit for the different data sets we consider for the DY 3 model, as calculated by smelli-2.2.0 for MX = 3 TeV. We display the total χ 2 for each data set along with the number of observables n and the data set's p−value. The set named 'quarks' contains BR(Bs → φµ + µ − ), BR(Bs → µ + µ − ), ∆ms and various differential distributions in B → K ( * ) µ + µ − decays among others, whereas 'LFU FC-NCs' contains R K ( * ) and some B meson decay branching ratios into di-taus. Our sets are identical to those defined by smelli-2.2.0 and we refer the curious reader to its manual [67], where the observables are enumerated. We have updated RK and BR(B s,d → µ + µ − ) with the latest LHCb measurements as detailed in §4. Two free parameters of the model were fitted: θ23 = −0.181 and gX = 0.253.
vide the driving force and so we neglect to show them 12 .
We shall now turn to the DY 3 model fit results, where these comments about flavour observables also apply.

DY 3 model fit results
We summarise the quality of the fit for the DY 3 model at the best-fit point, for M X = 3 TeV, in Table 7. We  TeV. Shaded regions are those preferred by the data set in the legend at the 95% confidence level (CL). The global fit is shown by the solid curves, where the inner (outer) curves show the 70% (95%) CL regions, respectively. The set named 'quarks' contains BR(Bs → φµ + µ − ), BR(Bs → µ + µ − ), ∆ms and various differential distributions in B → K ( * ) µ + µ − decays among others, whereas 'LFU FCNCs' contains R K ( * ) and some B meson decay branching ratios into di-taus. Our sets are identical to those defined by smelli-2.2.0 and we refer the curious reader to its manual [67], where the observables are enumerated. We have updated RK and BR(B s,d → µ + µ − ) with the latest LHCb measurements as detailed in §4. The black dot marks the locus of the best-fit point. see a much improved fit as compared to the SM (by a ∆χ 2 = 39) and a similar (but slightly worse) quality of fit compared to the Y 3 model, as a comparison with Table 6 shows.
The constraints upon the parameters θ 23 and g X are shown in Fig. 11. Although the figure is for M X = 3 TeV, the picture remains approximately the same for 2 < M X /TeV < 10. We see that, as is the case for the Y 3 model, there is a region of overlap of the 95% CL regions of all of the constraints.
The pulls in the EWPOs for the best-fit point of the DY 3 model are shown in Fig. 12. Like the Y 3 model above, we see a fit comparable in quality to that of the SM. Again, the DY 3 model predicts M W to be a little higher than in the SM, agreeing slightly better with the experimental measurement.
The behaviour of the fit in various directions in parameter space around the best-fit point is shown in Figs. 13-15. Qualitatively, this behaviour is similar to that of the Y 3 model: the EWPOs and NCBAs imply that g X should not become too large. The mixing observable ∆m s prevents −θ 23 from becoming too large, and the fits are insensitive to M X varied between 2 TeV and 10 TeV so long as g X is scaled linearly with M X . Fig. 16 shows various pulls of interest at the best-fit point of the DY 3 model. Better fits (than the SM) to sev-      Table 8. Goodness of fit for the different data sets we consider for the original DY3 model, as calculated by smelli-2.2.0 for MX = 3 TeV. We display the total χ 2 for each data set along with the number of observables n and the data set's p−value. The set named 'quarks' contains BR(Bs → φµ + µ − ), BR(Bs → µ + µ − ), ∆ms and various differential distributions in B → K ( * ) µ + µ − decays among others, whereas 'LFU FC-NCs' contains R K ( * ) and some B meson decay branching ratios into di-taus. Our sets are identical to those defined by smelli-2.2.0 and we refer the curious reader to its manual [67], where the observables are enumerated. We have updated RK and BR(B s,d → µ + µ − ) with the latest LHCb measurements as detailed in §4. Two free parameters of the model were fitted: θ23 = 0.122 and gX = 0.428.
eral NBCA observables are partially counteracted by a worse fit to the ∆m s observable.

Original DY 3 model fit results
We display the overall fit quality of the original DY 3 model in Table 8. By comparison with Tables 1 and 6 we see that although its predictions still fit the data significantly better than the SM (∆χ 2 is 32), the original DY 3 model does not achieve as good fits as the other models. For the sake of brevity, we have refrained from including plots for it 13 . Instead, it is more enlightening to understand the reason behind this slightly worse fit, which is roughly as follows. The coupling of the X boson to muons in the original DY 3 model is close to vector-like, viz. L = g X /6(µ / X(5P L + 4P R ))µ + . . . (where P L , P R are left-handed and right-handed projection operators, respectively), which is slightly less preferred by the smelli-2.2.0 fits than an X boson coupled more strongly to left-handed muons [67]. This preference is in large part due to the experimentally measured value of BR(B s → µ + µ − ), which is somewhat lower than the SM prediction [4][5][6][7], and is sensitive only to the axial component of the coupling to muons. Compared to the Y 3 model and the DY 3 model, the fit to BR(B s → µ + µ − ) is worse when the DY 3 model fits other observables well. The p−value is significantly lower than the canonical lower bound of 0.05, indicating a somewhat poor fit. well as non-zero C 2233 eu in the case of the DY 3 model and the DY 3 model) can give a reasonable fit to the NCBA data [89] via a W -boson loop as in Fig. 17 14 Such a scenario would require that V ts originates from mixing entirely within the up-quark sector. This qualitatively different quark mixing ansatz therefore provides a motivation to consider the θ 23 = 0 scenario separately. In the θ 23 = 0 limit, we have that Λ  Tables 4 and 5. We note that contributions to the NCBAs arising from W loops such as the one in Fig. 17 are nevertheless always included (through renormalisation group running) by smelli-2.2.0, even for θ 23 = 0.
While it is true that much of the tension with the NCBAs can be ameliorated by such W loop contributions, we find from our global fits that the corresponding values for g 2 X /M 2 X are far too large to simultaneously give a good fit to the EWPOs in this θ 23 = 0 limit. As our results with a floating θ 23 suggest (see e.g. Fig. 6 along θ 23 = 0), as far as the EWPOs go, SM-like scenarios are strongly preferred, since EWPOs quickly exclude any region that might resolve the tension with NCBAs. This is even more so than for the simplified model studied in [89], where already a significant tension with Z → µµ was pointed out. In our case, besides several stringent bounds from other observables measured at the Z 0 -pole for g X ≈ 1 and M X ≈ O(3) TeV (which predicts C 2223 lq just large enough to give a slightly better fit to R K and R K * than the SM) the predicted W −boson mass is more than 5σ away from its measured value. This stands in stark contrast with the overall better-than-SM fit we find for M W when θ 23 = 0.
14 The connection between the EWPOs and minimal flavourviolating one-loop induced NCBAs has recently been addressed in Ref. [90], encoding the inputs in the SMEFT framework. Third family hypercharge models are in a sense orthogonal to this analysis, since they favour particular directions in nonminimal flavour violating SMEFT operator space generated already at the tree-level. The possibility of accommodating the NCBAs through such loop contributions was first pointed out in Ref. [91] with a more general study for b → s transitions presented in Ref. [92].

Discussion
Previous explorations of the parameter spaces of third family hypercharge models [41,43] capable of explaining the neutral current B−anomalies showed the 95% confidence level exclusion regions from various important constraints, but these analyses did not include electroweak precision observables. Collectively, the electroweak precision observables were potentially a model-killing constraint because, through the Z 0 − Z mixing predicted in the models, the prediction of M W in terms of M Z is significantly altered from the SM prediction. This was noticed in Ref. [66], where rough estimates of the absolute sizes of deviations were made. However, the severity of this constraint on the original Third Family Hypercharge (Y 3 ) Model parameter space was found to depend greatly upon which estimate 15 of the constraint was used. In the present paper, we use the smelli-2.2.0 [67] computer program to robustly and accurately predict the electroweak precision observables and provide a comparison with empirical measurements. We have thence carried out global fits of third family hypercharge models to data pertinent to the neutral current B−anomalies as well as the electroweak precision observables. This is more sophisticated than the previous efforts because it allows tensions between 217 different observables to be traded off against one another in a statistically sound way. In fact, at the best-fit points of the third family hypercharge models, M W fits better than in the SM, whose prediction is some 2σ too low.
One ingredient of our fits was the assumption of a fermion mixing ansatz. The precise details of fermion mixing are expected to be fixed in third family hypercharge models by a more complete ultra-violet model. This could lead to suppressed non-renormalisable operators in the third family hypercharge model effective field theory, for example which, when the flavon acquires its VEV, lead to small mixing effects. Such detailed model building seems premature in the absence of additional information coming from the direct observation of a flavour-violating Z , or indeed independent precise confirmation of NCBAs from the Belle II [93] experiment. Reining in any urge to delve into the underlying model building, we prefer simply to assume a non-trivial structure in the fermion mixing matrix which changes the observables we consider most: those involving the left-handed down quarks. Since the neutral current B−anomalies are most sensitive to the mixing angle between left-handed bottom and strange quarks, we have allowed this angle to float. But the other mixing angles and complex phase in the matrix have been set to some (roughly mandated but ultimately arbitrary) values equal to those in the CKM matrix. We have checked that changing these arbitrary values somewhat (e.g. setting them to zero) does not change the fit qualitatively: a change in χ 2 of up to 2 units was observed. It is clear that a more thorough investigation of such variations may be- come interesting in the future, particularly if the NCBAs strengthen because of new measurements. We summarise the punch line of the global fits in Table 9. We see that, while the SM suffers from a poor fit to the combined data set, the various third family hypercharge models fare considerably better. The model with the best fit is the original Third Family Hypercharge Model (Y 3 ). We have presented the constraints upon the parameter spaces of the Y 3 model and the DY 3 model in detail in §4. The qualitative behaviour of the Y 3 model and the DY 3 model in the global fit is similar, although the regions of preferred parameters are different.
It is well known that ∆m s provides a strong constraint on models which fit the NCBAs and ours are no exception: in fact, we see in Figs 10,16 that this variable has a pull of 2.7σ (2.1σ) in the Y 3 model (DY 3 model), whereas the SM pull is only 1.1σ, according to the smelli-2.2.0 calculation. The dominant beyond the SM contribution to ∆m s from our models is proportional to the Z coupling tosb quarks squared, i.e. [g X (Λ ) 23 is adjustable because θ 23 is allowed to vary over the fit, and ∆m s provides an upper bound upon |g X (Λ (d L ) ξ ) 23 |. On the other hand, in order to produce a large enough effect in the lepton flavour universality violating observables to fit data, the product of the Z couplings tosb quarks and to µ + µ − must be at least a certain size. Thus, models where the Z couples more strongly to muons because their U (1) X charges are larger fare better when fitting the combination of the LFU FCNCs and ∆m s . The Z coupling to muons is 1/2 for the Y 3 model and 2/3 for the DY 3 model, favouring the DY 3 model in this regard.
Despite the somewhat worse fit to ∆m s for the Y 3 model as compared to the DY 3 model, Table 9 shows that, overall, the Y 3 model is a better fit. Looking at the flavour observables in detail, it is hard to divine a single cause for this: it appears to be the accumulated effect of many flavour observables in tandem. The difference in χ 2 between the Y 3 model and the DY 3 model of 4.4 is not large and might merely be the result of statistical fluctuations in the 219 data; indeed 1.2 of this comes from the difference of quality of fit to the EWPOs.
All of the usual caveats levelled at interpreting p−values apply. In particular, p−values change depending upon exactly which observables are included or excluded. We have stuck to pre-defined sets of observables in smelli-2.2.0 in an attempt to reduce bias. However, we note that there are other data that are in tension with SM predictions which we have not included, namely the anomalous magnetic moment of the muon (g − 2) µ and charged current B−anomalies. If we were to include these observables, the p−values of all models in Table 9 would lower. Since third family hypercharge models give the same prediction for these observables as the SM, each model would receive the same χ 2 increase as well as the increase in the number of fitted data. However, since our models have essentially nothing to add to these observables compared with the SM, we feel justified in leaving them out from the beginning. We could have excluded some of the observables that smelli-2.2.0 includes in our data sets (obvious choices include those that do not involve bottom quarks, e.g. K ) further changing our calculation of the p−values of the various models.
As noted above, as far as the third family hypercharge models currently stand, the Z contribution to (g − 2) µ is small [28]. In order to explain an inferred beyond the SM contribution to (g − 2) µ compatible with current measurements ∆(g − 2) µ /2 ≈ 28 ± 8 × 10 −10 , one may simply add a heavy (TeV-scale) vector-like lepton representation that couples to the muon and the Z at one vertex. In that case, a one-loop diagram with the heavy leptons and Z running in the loop is sufficient and is simultaneously compatible with the neutral current B−anomalies and measurements of (g − 2) µ [28].
Independent corroboration from other experiments and future B−anomaly measurements are eagerly awaited and, depending upon them, a re-visiting of global fits to flavour and electroweak data may well become desirable. We also note that, since electroweak precision observables play a key rôle in our fits, an increase in precision upon them resulting from LHC or future e + e − collider data, could also prove to be of great utility in testing third family hypercharge models indirectly. Direct production of the predicted Z [43,65] (and a measurement of its couplings) would, along with an observation of flavonstrahlung [40], ultimately provide a 'smoking gun' signal.