Third-Family Quark-Lepton Unification and Electroweak Precision Tests

We analyze the compatibility of the hypothesis of third-family quark-lepton unification at the TeV scale with electroweak precision data, lepton flavor universality tests, and high-$p_T$ constraints. We work within the framework of the UV complete flavor non-universal 4321 gauge model, which is matched at one loop to the Standard Model Effective Field Theory. For consistency, all electroweak precision observables are also computed at one loop within the effective field theory. At tree level, the most sizeable corrections are to $W\rightarrow \tau\nu_\tau$ and $Z \to \nu_\tau \nu_\tau$ due to integrating out a pseudo-Dirac singlet fermion required by the model for neutrino mass generation. At loop level, the new colored states of the model generate large flavor-universal contributions to the electroweak precision observables via leading- and next-to-leading log running effects, yielding a significant improvement in the electroweak fit (including an increase in the $W$-boson mass). These effects cannot be decoupled if the model addresses the charged-current $B$-meson anomalies. Overall, we find good compatibility between the data sets, while simultaneously satisfying all low- and high-energy constraints.

of lepton flavor universality (LFU) violations in the b → slℓ [23] system changes little concerning the theoretical motivation of and interest in the model, although there is clearly less experimental evidence for a large NP effect in light-family leptons.
The rich phenomenology following from the quark-lepton unification hypothesis has been studied in detail in several previous works, including implications for τ decays [24], ∆F = 2 and ∆F = 1 transitions [25][26][27], and high-p T collider signatures [28][29][30][31]. An element that has been missing so far, however, is a systematic study of the implications for electroweak precision observables (EWPO), which constitutes the main objective of this paper.
The interplay between semi-leptonic interactions at the TeV scale, generated by flavor non-universal interactions, and EWPO has already been pointed out in a pure Effective Field Theory (EFT) context [32,33]: the running of the (non-universal) four-fermion operators into Higgs-fermion operators results in sizable corrections to the EWPO. Having a UV-complete model at hand, we now go beyond these initial studies including all the relevant finite pieces resulting from the matching of the 4321 model into the EFT of the SM (the so-called SMEFT). In this respect, we analyze for the first time in a systematic way the effects of new fermion fields (and related interactions), which are present in all realistic scenarios featuring third-family quark-lepton unification at the TeV scale. First, these models necessarily include vector-like fermions in order to describe the mixing between the light and heavy SM families. Second, an inevitable consequence of unifying quarks and leptons of the third family is the prediction of a degenerate top quark and tau neutrino. In order to resolve this issue, an inverse-seesaw mechanism with additional gauge singlet fermions is a necessary ingredient [2]. In addition to the U 1 LQ, the model requires the existence of two massive neutral vector bosons. As we shall show, all of these new states generate important effects in EWPO. For example, the exotic fermions in many cases induce one-loop matching contributions that can compete in size with pure renormalization group (RG) effects. Additionally, four-quark operators generated after integrating out the neutral vectors can run beyond the leading-log approximation into relevant contributions for the electroweak (EW) fit. Besides providing analytical expressions to evaluate all these effects in general terms, we perform a detailed phenomenological analysis of different EWPO, identifying the most interesting ones to further test this hypothesis in the future. Finally, we analyze the compatibility of EWPO with LFU tests in b → cτ ν and τ -decays, as well as with high-p T constraints by performing a global fit to the data.
The paper is organized as follows. In Section II we define the non-universal 4321 model. In Section III we discuss the matching procedure from the non-universal 4321 model to the SMEFT, the main running effects from the 4321 breaking scale down to the electroweak scale, and the one-loop computation of the EWPO in the EFT. The phenomenological discussion as well as the results of the global fit are presented in Section IV and summarized in the Conclusions. A series of appendices contain more technical aspects: a detailed discussion about one-loop matching contributions (A); a discussion and analytical expressions for the RG running of the Wilson coefficients (WCs) (B); one-loop expressions for the EWPO in terms of SMEFT coefficients (C); phenomenological expressions for the LFU ratios in b → cℓν transitions (D).

II. THE MODEL
As already anticipated, the gauge group of the model is where SU (3) c is identified as the diagonal subgroup of SU (4) h × SU (3) l and the flavor-universal group SU (2) L acts like in the SM. The hypercharge is given by Y = X + 2/3 T 15 4 , where T 15 4 = 1 2 √ 6 diag(1, 1, 1, −3) is a generator of SU (4) h , meaning that U (1) X coincides with hypercharge for the light families (see Table I for the full matter field content). Without mass mixing, the three SU (4) h multiplets ψ L and ψ ± R can be identified with the SM thirdgeneration fermions (with the addition of a right-handed neutrino). The only exotic fermions are a vector-like fermion χ L,R = (Q L,R , L L,R ), with the same gauge quantum numbers as ψ L , and a gauge singlet S L needed to achieve an acceptable tau neutrino mass.
The scalar fields Ω 1 , Ω 3 , and Ω 15 mediate the G 4321 → SM breaking at the TeV scale. 1 This results in three heavy gauge fields, transforming under the SM as U 1 ∼ (3, 1, 2/3), G ′ ∼ (8, 1, 0) and Z ′ ∼ (1, 1, 0). More details about Field SU (4) h SU (3) l SU (2)L U (1)X symmetry breaking can be found in [26]. In the G 4321 broken phase, a mixing between chiral and vector-like fermions is generated. Defining Ψ ′⊺ q = (q ′1 L q ′2 L q ′3 L Q ′ L ) and Ψ ′⊺ ℓ = (ℓ ′1 L ℓ ′2 L ℓ ′3 L L ′ L ), the mass mixing can be written as where, without loss of generality, the mass vector can be decomposed as with m Q and m L being the physical vector-like fermion masses. Here, O q,ℓ describe the mixing between different SU (4) h representations (q ′2 L − Q ′ L and ℓ ′2 L − L ′ L ), parameterized by the angle θ q,ℓ , while W q,ℓ parameterizes the mixing amongst SU (4) h states (q ′3 L − Q ′ L and ℓ ′3 L − L ′ L ). After moving to the mass basis, Ψ ′ q,ℓ → W † q,ℓ O † q,ℓ Ψ q,ℓ , the relevant interactions of the fermions with the massive gauge bosons read where g 4 indicates the SU (4) gauge coupling, ϕ R is the relative phase between left-and right-handed U 1 currents, and G ′ µ = G ′a µ T a with T a being the Gell-Mann matrices. The flavor structure is summarised in the matrices where we have neglected terms O(g 2 s /g 2 4 ) and O(g 2 Y /g 2 4 ), where g s and g Y are the SM strong and hypercharge couplings. By re-phasing q 3 L , ℓ 3 L , and the left-handed vector-like fermion fields, the matrix W = W q W † ℓ can be written as a 2 × 2 real rotation, parameterized by an angle χ, Similarly, we can write the Yukawa interactions in the mass basis as where the Yukawa vectors can be expressed as with y t,b being the top and bottom Yukawa couplings, respectively, and Y ± = |Y ± |e iϕ± a complex parameter. Stringent constraints from B s −B s mixing mediated by the neutral G ′ , Z ′ gauge bosons at tree-level imply that a significant amount of alignment to the down-quark mass basis in the 2-3 sector is phenomenologically required [27]. The dominant breaking of U (2) q that generates 2-3 mixing in the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix must therefore be realized in the up-quark sector, corresponding to the limit Y − ≪ Y + . In this case, the contribution of Y − to the CKM matrix element V cb can be neglected, and we have Working in the limit Y − → 0, the leading Yukawa interactions with the Higgs field are those involving y t and Y + where the couplings y ν and Y ν + are defined as 2 Finally, the smallness of the tau-neutrino mass is obtained via the addition of a singlet fermion S L in order to implement the inverse-seesaw mechanism [2,3]. Effectively, it interacts with ν 3 R through the Lagrangian resulting in a pseudo-Dirac pair split by the lepton number violating parameter µ. In particular, in the limit µ ≪ λ R ⟨Ω 1 ⟩, these states have nearly degenerate masses m R ≈ λ R ⟨Ω 1 ⟩ ± µ/2 tied to the 4321 breaking scale. After EW symmetry breaking there is also a light, active Majorana state with mass m ν ≈ µ(y ν ⟨H⟩/m R ) 2 . In order to obtain small active neutrino masses for m R ∼ TeV and y ν ∼ O(1), we need to work near the (technically natural) U (1) L preserving limit µ ∼ 0. Thus, it is a very good approximation to consider only a single heavy Dirac neutrino of mass m R , and we work in this limit in what follows.

III. MATCHING, RUNNING, AND OBSERVABLES IN THE SMEFT
In this section, we discuss the matching of the model introduced in the previous section to the SMEFT. We normalize the Lagrangian as where the operators of interest, O k , are listed in Table II and III. We perform tree-level and one-loop matching, with the following guiding principles: 2 In the model as written, the W -matrix is the only source of quark-lepton splitting and the angle χ and the phase ϕ R = −arg(yτ /y b ) are therefore fixed by yτ /y b . However, on general grounds, we expect additional sub-leading SU (4) h -breaking corrections such as ψ L Ω 15 Hψ − R from integrating-out heavier fields that do not impact the low-energy phenomenology. These corrections can however have a significant impact on the (small) down-type Yukawas, so we treat χ and ϕ R as free parameters in what follows.  • Tree-level. All operators generated by tree-level integration of the heavy fields are included, except those that do not enter EWPO, high-p T constraints, or b → cτ ν observables, neither directly nor through RGE effects. In particular, integrating out the vector-like quark Q and the right-handed neutrino ν R yields some of the Higgsfermion current operators in Table II, while integrating out the heavy gauge bosons (U 1 , G ′ , and Z ′ ) gives rise to the four-fermion operators in Table III.
• One-loop. One-loop effects are accounted for by taking corrections from y t , Y + , and g 4 for the operators entering the electroweak fit directly. In practice, this means that we give the one-loop matching conditions for C (1,3) Hℓ , C (1,3) Hq , and C HD , with third-generation flavor indices only. We do not include one-loop matching for [O Hu ] 33 since it impacts EWPO only via RG mixing into O HD . Operators with light fermion fields are always suppressed, e.g. by s q,ℓ ≲ 0.1, and have a minimal impact on the electroweak fit. Finally, we include one-loop corrections in g 4 for the semi-leptonic operators relevant for b → cτ ν and high-p T observables.
In the following, we give the results for the tree-level matching and the leading running effects. The expressions for the finite parts from the one-loop computation and a more detailed discussion of the matching procedure are given in Appendix A.
A. Tree-level matching All tree-level matching conditions are given in Table IV. For example, integrating out the vector-like quark Q gives [O Hu ] 33 , whose RG mixing into the O HD operator (proportional to y 2 t ) provides a shift to the W mass and couplings of the EW gauge bosons to fermions. Furthermore, integrating out the massive Dirac neutrino gives [C Hℓ ] 33 , resulting in tree-level modifications to Z-boson couplings to τ -neutrinos as well as the W -boson coupling to τ ν τ . In the matching of the four-fermion operators, we introduce the effective scale while the coupling matrix β L (omitting couplings with the vector-like fermions) explicitly reads where the fact that the U 1 does not interact with the first generation prior to electroweak symmetry breaking is manifest. Of the four-fermion operators generated by integrating out the G ′ and Z ′ , we retain only those with purely and the Wilson coefficient [C Hℓ ] 22 enters the electroweak fit via a modification of the muon decay width that is used to extract the Fermi constant G F . However, including s ℓ as a parameter would open up a new flat direction in the EW fit, requiring an s ℓ -sensitive observable to lift. In this work, we set s ℓ = 0 due to a lack of experimental evidence for large LFU violating effects in the b → slℓ system.

B. Leading running effects
Working with the complete UV model, we are able to compute one-loop corrections arising from the matching conditions at the scale of new physics. The leading effects are encoded in the leading-log (LL) and next-to-leading log (NLL) contributions arising due to the RG evolution of those operators induced at the tree level. In this work, we consider LL running in the three largest SM couplings y t , g s , g L , as well as NLL running in y t and g s . Here, we summarize the leading running effects, which are due to y t .

Leading log
In the following, we report the LL contributions to the Wilson coefficients of the operators affecting EWPO. Starting with O qq , and Cuu, as defined in Table IV.
By SU (2) L arguments, the contribution to the Wilson coefficient of the triplet operator can be obtained as Moreover, the running of the semi-leptonic operators O Hℓ , as already shown in [24]. In addition, the operators O (1,3) Hℓ also run into themselves. The complete leading-log result is therefore and [C 2. Next-to-leading log We find that the NLL running can give important contributions to the operators relevant for the electroweak fit. It can be calculated as where [C] represents the vector of all dimension-6 Wilson coefficients and A is the anomalous dimension matrix that can be read from [34][35][36] (see Appendix B for more details). The most important contribution comes from the top Yukawa-induced NLL running of the four-quark operators O qq , and O uu generated by tree-level coloron exchange into O HD , corresponding to the two-loop diagram shown on the right in Fig. 1. This NLL contribution to C HD reads where the WCs are understood to have purely third-family flavor indices. To illustrate the potential importance of this effect, we convert the leading and next-to-leading log contributions to C HD into shifts in m W (Eq. (C22)) and plot the result in Fig. 2. One can see for Λ U ≲ 2 TeV that the NLL running can easily be of the same order or larger than the LL contribution. Additionally, both the LL and NLL contributions to C HD have the same sign and thus add coherently. We include this contribution as well as all NLL running contributions in the couplings g s and y t from third-family operators generated at the tree level in the UV in our analysis. The complete expressions for the NLL effects are given in Appendix B 2, and the proper treatment of the top Yukawa y t is discussed in detail in Appendix B.
As a final comment, we note that O HD violates custodial symmetry, and therefore its generation should correspond to custodial violations within the model. In the case of the LL contribution, the violation comes from Y − ≪ Y + , which is required phenomenologically to pass the bounds on B s −B s mixing. In the case of the NLL contribution, the custodial symmetry violation is of SM origin, namely y b ≪ y t . It is interesting to note that taking the symmetric limit in both cases is not possible, so these effects cannot be switched off by invoking custodial symmetry. We are varying the relevant parameters in the ranges mQ ∈ [1.5, 3] TeV and m G ′ ∈ [3, 3.5] TeV. In the darker blue regions, the IR scale is fixed to µEW = mt. In the right plot, the light-blue region shows the impact of varying the IR scale from mt/2 to 2mt. The gray horizontal line corresponds to the ∆mW necessary to accommodate the averaged experimental value mW = 80.379 GeV before CDF II (see Eq. (38)).

C. Observables at one-loop accuracy
Besides the tree-level matching and running effects, we include the one-loop matching in y t , Y + , and g 4 of the operators involved in the EW fit, as given in Appendix A. Once this is done, for consistency we have to include the one-loop corrections in y t of the EWPO. These corrections are of the same size as the non-log-enhanced terms of the one-loop matching. If we include these finite pieces from the one-loop matching, the computation of the EWPO at the treelevel in SMEFT would be basis dependent due to the existence of evanescent operators [37]. However, the one-loop calculation of the observables cancels these ambiguities. Moreover, the IR scale dependence of the observables is also canceled at the leading-log order. For all beyond leading-log effects, we choose to fix µ EW = m t . Working in the {α EM , m Z , G F } input scheme [38], the EWPO can be written in terms of the W and Z bosons vertex corrections and W boson mass shift. In Appendix C we provide all the expressions for these pseudo-observables, including one-loop corrections in y t .

IV. PHENOMENOLOGY AND GLOBAL FIT
In the following, we discuss the phenomenology of the model given in Section II featuring third-family quark-lepton unification. We begin with some general remarks about the expected NP effects in EWPO, breaking them down into flavor non-universal and flavor universal effects: • Flavor non-universal effects. At tree level, the only operators generated by the model that directly affect EWPO are [C (1) Hℓ ] 33 and [C Hℓ ] 33 , leading to effects in EW gauge boson couplings to third-family leptons. However, because they are generated with the relation C (1) Hℓ , there is no modification to Z-boson couplings to tau leptons at tree level. Indeed, only W → τ ν τ and Z → ν τ ν τ are affected at tree level, as can be seen in Appendix C. The sign of the effect is such that existing tensions in W → τ ν τ are worsened, while Z → ν τ ν τ is made more compatible with data. At the loop level, this statement remains true when considering only LL running in y t (see Eq. (21)), so the leading breaking to the relation C (1) Hℓ comes from LL running effects in the SU (2) L gauge coupling g L . We give the expressions for this effect in Appendix B 3 and take it into account in our analysis. However, the effect in Z → τ L τ L is small and we find no significant constraint from it. Similarly, no modification of Z-boson couplings to bottom quarks is generated at tree level. The leading effect in Z → b L b L comes from the breaking of the relation C (1) Hq by the y t -induced running of C Hu into C (1) Hq (Eq. (19)), LL running in g L , and NLL running in y t and g s . These contributions are all of a similar size and we take them into account. Again, we find no significant correction to Z → b L b L vertex. There is also a coloron-induced NLL running effect leading to a sizeable modification of the Z-boson coupling to right-handed bottom quarks (B11), but this coupling is poorly constrained and leads to no significant impact on the fit.
• Flavor universal effects. As discussed in detail in Section III B, LL and NLL running of tree-level operators induced by the new colored states of the model lead to large flavor universal effects in C HD (see Fig. 1). Furthermore, these effects have the same sign and thus add coherently, giving a large universal NP effect that is preferred by the current EW data. In particular, the sign of the effect is such that it increases m W as well as Z-pole observables such as Γ Z , the effective number of neutrinos N eff ν , and asymmetry observables such as A ℓ,q and A FB b,c via flavor universal shifts δ U (T 3 , Q, C HD ) to the Z-boson couplings (shown in Appendix C). Finally, the flavor non-universal operator [O Hq ] ii via LL runnning in g L . While this effect is irrelevant compared to the universal contribution from C HD , we take it into account via the expressions in Appendix B 3.
As we will see below, this flavor universal effect in C HD , together with the tree-level non-universal effects in W → τ ν τ and Z → ν τ ν τ , are dominantly responsible for the behavior of the EW fit.
We now proceed with a global fit taking into account data from low-energy experiments, EWPO, and neutral current Drell-Yan processes at high-p T . In principle, the parameters relevant for the phenomenological discussion are: and m Z ′ are the masses of the heavy 4321 gauge bosons U 1 , G ′ , and Z ′ , respectively. Notice that s q can be written in terms of Y + using Eq. (11), which eliminates it as a free parameter. While Y + = |Y + |e iϕ+ is in general a complex parameter, a non-vanishing phase would make it necessary to re-phase the SM fields to match the standard CKM phase convention (see Appendix B of [27] for details). This re-phasing would affect the couplings of the U 1 LQ to second-family quarks, namely (β L ) 2i → e −iϕ+ (β L ) 2i in Eq. (17). To obtain the correct sign and maximize the effect of the LQ on b → cτ ν observables, we fix ϕ + = 0 and consider Y + to be real in what follows. Furthermore, we fix s ℓ = 0 according to the discussion below Eq. (18). As can be seen in Table IV, all tree-level effects are described by the effective scales up to mass splittings between the 4321 gauge bosons encoded by x Z ′ , x G ′ ∼ 1. This is no longer true beyond tree level, where the masses appear alone inside logs and loop functions. Still, the leading behavior is captured by the ratios in Eq. (25), so our strategy will be to fix the masses m U , x G ′ , x Z ′ , m L , and m R to natural values while allowing Λ U , Y + , and m Q to vary in the fit freely. The reason we do not fix all masses is that, from Eq. (13), the coupling is not an independent parameter. Furthermore, we choose to vary m Q instead of m R since the scale Λ R is bounded from below by the process W → τ ν in the EWPO and τ -decays, so fixing m R prevents the fit from finding points with very light m R by fine-tuning the two terms in Eq. (26). The benchmark point for the fixed parameters, and the variables freely varied in the fit, are summarized as follows

A. Fit Observables and Methodology
We construct the likelihoods used in the fit by building the χ 2 -function, defined as where σ −2 is the inverse of the covariance matrix. Here, we discuss all observables included in the fit and the corresponding theory predictions. The observables we choose to include are: • LFU tests in b → cτ ν transitions. These include the ratios R D and R D * , as well as R Λc , defined as The expressions for shifts of these observables from their SM predictions in terms of SMEFT operators are given in Appendix D. The world averages of R D ( * ) ratios, including the recent LHCb measurement, read [39] R exp with correlation ρ = −0.29. For the SM predictions, we use the HFLAV averages [39] R SM Additional discussion about the SM predictions for R D ( * ) and their uncertainties can be found in [40][41][42][43][44][45].
Regarding R Λc , we use the recent LHCb measurement as the input [46] and the result in [47] as the SM prediction • LFU tests in τ decays. These are expressed through the ratios expected to be equal to one in the SM by definition. We use the experimental average from [48] g τ g e,µ ℓ,π,K = 1.0012 ± 0.0012 .
The theory predictions for the τ LFU tests can be found in Appendix D.
• Electroweak precision observables. At the Z-and W -pole, we include all the observables listed in Tables 1  and 2 of Ref. [38] using the same SM theory predictions, and we follow the same input scheme. For the NP theory predictions, we parameterize deviations from the SM due to NP effects via shifts in the W and Z boson couplings, δg W,Z k , and the shift in the W -boson mass δm W . To simplify notation, we include δm W in the vector δg W,Z k and define the theory prediction for a given EWPO as The theory predictions for the δg Z,W k in terms of the SMEFT Wilson coefficients at one-loop in y t , as well as the map between the EWPO and the δg Z,W k that determines the α ik , are given in Appendix C.
We will discuss and compare the fit in the two scenarios, one with and one without the latest CDF II result on m W [49]. In the former case, given the statistical incompatibility of the CDF II measurement with all the previous determinations, a strategy needs to be defined in order to be able to perform a combination. Our choice is to assume that all experimental errors have been underestimated. Penalizing all measurements "democratically", and inflating the errors until the χ 2 per degree of freedom is equal to one, one finds m W = 80.410 ± 0.015 GeV for the averaged W mass. Therefore, we consider two scenarios with different input values for m W m old W = 80.379 ± 0.012 GeV without CDF II , • High-p T constraints. We consider the tail of pp → τ τ distributions obtained at the LHC experiments. The likelihood is obtained using HighPT [50], and for the theory input we use the semi-leptonic SMEFT coefficients computed at NLO in g 4 , as discussed in Appendix A 6. The likelihood is constructed using data from ATLAS [51]. However, it is interesting to consider also the recent search by CMS [52] which shows a small excess, leading to a weaker bound on the NP scale. As done in Ref. [53], we rescale the likelihood obtained from HighPT to match the NLO predictions derived in Ref. [54] for the ATLAS and CMS searches in the b-tag channel. In addition, a lower bound on the vector-like quark mass m Q is implemented by introducing a term in the likelihood of the form where m bound Q = 1.5 TeV is the lower bound on the mass at 95% CL [55] from direct searches for SU (2) L -doublet vector-like quarks pair-produced via QCD and decaying dominantly to Ht and Zt, as we expect in this model.

B. Fit Results
As discussed in Section IV A, we have the likelihoods χ 2 b→cτ ν , χ 2 EWPO , χ 2 τ -LFU , and χ 2 high-p T . We define the global likelihood simply as the sum of the individual likelihoods Fixing the parameters as discussed in the beginning of Section IV, the global likelihood is a function of 5-parameters, namely χ 2 (Λ U , Y + , m Q , χ, ϕ R ). However, as can be seen in Fig. 3, when choosing m new W the global χ 2 -function is very flat for χ ∈ [35 • , 80 • ], where it changes by less than 2 units along the red line giving the optimal value of Λ U . However, the optimal value of Y + (dashed black contours) is not flat in this range. Instead, it decreases from around 1 for χ = 35 • to around 0.4 for χ = 80 • . Generally speaking, larger values of Y + favor the EWPO which calls for non-zero C HD ∝ |Y + | 2 , while smaller values of Y + favor the b → cτ ν observables, which call for larger s q ∝ |Y + | −1 . Due to the flatness of the χ 2 -function between 35 • and 80 • , we choose to fix χ = 60 • for everything that follows, to allow for an intermediate value of Y + ∼ 0.6 that strikes a compromise between EWPO and b → cτ ν observables while making essentially no difference in the goodness of the fit. We have explicitly checked that this choice for χ is also good for the fit using m old W . Having fixed the angle χ = 60 • , we minimize the global χ 2 -function with respect to the remaining 4 parameters and report their best-fit values and 1σ-confidence intervals in Table V. The best-fit point (BFP) corresponds to ∆χ 2 = χ 2 SM − χ 2 BFP = 12.3 (15.4) in the case of m old W (m new W ), indicating a significant improvement over the SM hypothesis. 3 Next, we show the 1σ and 2σ preferred regions in the Λ U vs. Y + plane in Fig. 4 for the b → cτ ν observables (blue), EWPO (red), as well as the global preferred regions including all observables (blue lines without filling). The region below the dashed black line is excluded by CMS pp → τ τ data, while the gray shaded region is excluded by τ -LFU tests, both at 95% CL. In addition, the explicit relation between the couplings Y + and y ν , fixed by  Table V.
Eq. (26), allows us to show the preferred regions for y ν (upper horizontal axis in Fig. 4) simultaneously. Importantly, the parameters that are not varied in the plots are fixed to their global best-fit points, except in the case of m old W where we fix m Q = 3.2 TeV. This is due to the insensitivity of the fit to m Q when m Q ≳ 2.3 TeV, so we simply choose a theoretically reasonable value in the m old W case. We now offer some general conclusions from the results of the fit. In the case of b → cτ ν transition observables, we remark that the preferred region allows for larger Λ U values as |Y + | gets smaller. This behavior is mainly dictated by the fact that LFU ratios R D ( * ) and R Λc receive contributions proportional to β 23 L = s q s χ , where s q = V cb y t /|Y + |. Smaller Y + , therefore, increases this contribution, allowing for a higher effective scale Λ U to achieve the same size NP effect. Moreover, as shown in [27], K-K mixing and measurements of B(K + → π + νν) limit the magnitude of s q ≲ 0.1, thus requiring larger values of |Y + | > ∼ 0.35. On the other hand, the 1σ confidence interval for Y + translates to s q ∈ [0.04, 0.08]. It is very interesting that the results of this fit, which does not include kaon observables, are well compatible with the bound on the leading U (2) q breaking parameter s q found in Ref. [27]. Finally, moving in the direction of larger |Y + | (smaller s q ), we observe the simultaneous decrease in the value of Λ U necessary to accommodate the b → cτ ν data as expected.
Regarding τ -LFU tests, they exclude |y ν | ≳ 0.5 at 95% CL for almost all values of Λ U , with the bound becoming somewhat stronger for very low values (around 500 GeV) of Λ U . This is due to constructive addition between the tree-level (∝ −|y ν | 2 /m 2 R ) and one-loop (∝ c 2 Hℓ (see also the discussion below). These statements hold where the mass of the right-handed neutrino has been fixed to m R = 1.5 TeV. At the best-fit point, we have y ν = −0.13, and the 1σ confidence interval is y ν ∈ [−0.29, 0.06]. This corresponds to mixing angles θ 2 τ ≡ |y ν | 2 v 2 /(2m 2 R ) in the range θ 2 τ ∈ [5 × 10 −5 , 10 −3 ], where v = 246 GeV. As a consistency check for the fit, these values are well below the bound obtained from EWPO in the literature of θ 2 τ < 5.3 × 10 −3 [56,57]. Switching on only C (1) Hℓ = θ 2 τ /(2v 2 ), we find a similar but somewhat stronger bound of θ 2 τ < 2.3 × 10 −3 at 95% CL when combining the likelihoods from τ -LFU tests and EWPO that neatly explains the preferred range on |y ν | from the fit. We note that smaller values of m R will lead to a tighter upper bound on |y ν |, while the bound gets relaxed in the large m R limit. However, since m R ∝ ⟨Ω 1 ⟩, it cannot be fully decoupled without simultaneously decoupling the U 1 LQ. Therefore, when the scale of the model is low enough to address the charged-current B-anomalies, it always predicts sizeable unitarity violations, θ τ , in the active neutrino mixing matrix unless y ν is tuned small, as first pointed out in [2]. Coming next to the EWPO, one can observe that: i) When m old W is used as an EWPO (left panel in Fig. 4), the fit generally prefers lower values of |Y + | and Λ U ≳ 1 TeV, i.e. those regions in which m W does not receive a large shift. This is a consequence of the fact that δm W ∼ |Y + | 2 , from the LL contribution to C HD , as can be seen in Eq. (22). Additionally, there is the NLL contribution given in Eq. (24) that is proportional to Λ −2 U and adds constructively with the LL contribution. It is interesting that the lower bound of Λ U ≳ 1 TeV, coming dominantly from the NLL contribution to C HD , is stronger than that of τ -LFU tests when χ is fixed to 60 • .
ii) When m new W is used as an EWPO (right panel in Fig. 4), the EW fit prefers to have larger |C HD |, reflected in the fact that larger values of |Y + | are now preferred by EW precision data for Λ U ≳ 1.5 TeV, where the NLL effect is less relevant (as shown in Fig. 2). Furthermore, we note the preferred region also stretches to smaller values of |Y + | and Λ U , |Y + | ∈ [0.2, 0.6] and Λ U < 1.5 TeV. This is due to the NLL contributions to C HD from the four-quark operators induced by the coloron exchange at tree-level, Eq. (24). This effect is inversely proportional to Λ 2 U , and adds constructively with the NLL contribution, allowing a larger W mass to be accommodated even for small |Y + |, provided Λ U is small enough.
Finally, we note that for fixed values of the angle χ, the high-p T constraints are sensitive only to the value of Λ U . The black dashed line gives the 95% CL exclusion limit from CMS, which is providing by far the most stringent lower bound on the NP scale Λ U ≳ 1.25 TeV. Similarly, the 95% CL exclusion limit from the ATLAS data is Λ U ≳ 1.6 TeV. However, due to the fact that ATLAS currently has an under-fluctuation in the data and is therefore setting more stringent limits than expected (and also does not yet have a dedicated non-resonant search in the pp → τ τ channel), we consider the CMS bound of Λ U ≳ 1.25 TeV to be the current hard lower bound on the scale. Therefore, in the global fit, we have used the CMS likelihood for χ 2 high-p T . Finally, it is worth mentioning that the right-handed couplings of the U 1 LQ to the SM fermions can be reduced in less minimal versions of the model, e.g. via mixing with extra vector-like fermions. In such a case, the bound from high-p T constraints would be reduced, decreasing the tension with the current ATLAS data [53].

C. Selected Observables
Having discussed the general behavior of the fit, it is interesting to have a look at the observables that dominantly drive the fit. We define the pull of a particular observable O computed in a theoretical model M to be where σ exp is the experimental uncertainty for a given measurement O exp . In Table VI, we show the highest pulling observables for the best-fit point (BFP) of the NP model as well as the same observables within the SM. Furthermore, in Fig. 5 we compare the NP model predictions using m new W (blue regions) with the current experimental determinations for selected individual observables (red regions). We also show the 1σ interval for the model parameter being plotted in green, so the overlap of the blue and green regions should be considered as reasonably obtainable within the model. We can make the following observations: • R D /R D * [Figs. 5(a) and 5(b)]. At the best-fit point, we have δR D = 0.18 and δR D * = 0.04, so the model gives a very good fit for R D , while lying slightly below the central experimental value (≈ 1.7σ) for R D * . The reason why a large contribution to R D can be obtained even for relatively high Λ U is due to the right-handed b-τ coupling of the U 1 LQ in this model. As the fit prefers ϕ R = π, a large constructive scalar contribution to R D is generated via [C ℓedq ] 3332 (see D2). In general, as argued before, smaller values of |Y + | or Λ U are required to increase the contribution to both ratios, which then increases the tension with EWPO or high-p T , respectively. give a larger NP contribution via LL running into C HD . In Fig. 5(d), which shows δm W as a function of Λ U , one can see the impact of NLL running into C HD from coloron exchange, giving a large effect in m W for small values of Λ U . The small Λ U < 1.25 TeV region is, however, strongly in tension with pp → τ τ data. We note that the central value of m new W lies within the overlap of the green and blue regions, and is therefore achievable within the model. opposite directions, with NP effects causing more tension with τ -LFU measurements. As large Y + corresponds to negative y ν , the asymmetric shape of the N eff prediction for negative values of y ν stems from universal contributions to the Z-boson couplings to neutrinos via C HD . Such universal shifts cancel out in the LFU ratios and are therefore not seen in Fig. 5(h), where the model always increases the tension with data.

V. CONCLUSIONS
We have studied for the first time the phenomenological implications of a UV complete model featuring thirdfamily quark-lepton unification for electroweak precision tests. This was done by matching the UV model (defined in Section II) to the SMEFT at one loop, paired with a one-loop computation of the EWPO within the SMEFT (which is necessary to consistently capture all finite parts). We find good agreement between the model predictions for the EWPO and the experimental results, with an overall improvement of the EW fit with respect to the SM.
At tree level, the most sizeable corrections come from integrating out a heavy pseudo-Dirac singlet necessary to account for acceptable neutrino masses within the model. The primary impact of this state on EWPO and τ -LFU tests is to increase the pre-existing tension in the W -boson coupling to third-family leptons. However, one of the most important results of our analysis is the demonstration that higher-order effects play a key role. They give rise to a qualitatively different behavior of the fit compared to the inclusion of tree-level effects only, and break flat directions in the parameter space. In this respect, many of the results we have derived, from the matching to the SMEFT to the computation of EWPO, have a range of applicability that goes well beyond the specific model analyzed here. They provide an important addition to any phenomenological analysis of a wide class of models featuring vector-like fermions, extended gauge groups, or TeV-scale implementations of the inverse-seesaw mechanism.
Concerning the specific framework we have analysed, we find that the new colored states of the model, while not affecting EWPO at the tree level, generate a large flavor universal contribution at the loop level. This loop contribution comes via LL and NLL running, as well as finite one-loop matching contributions to the SMEFT operator O HD . This effect dominates the overall improvement in the electroweak fit via an increase in the W -boson mass and universal shifts in several other Z-pole observables. Furthermore, we observe good compatibility in the interplay of EWPO with other data, such as LFU tests at low energies and high-p T bounds, where the global best-fit point of the model gives a significant improvement over the SM hypothesis. In the case where the model addresses the charged-current B-anomalies, our combined analysis suggests the existence of new states not far above the TeV scale, giving effects in EWPO that cannot be decoupled. This sets a clear target for both current and near-future experiments, both at the intensity and high-energy frontiers. The one-loop diagrams in the full theory proceed through box diagrams propagating fermions and either vector or scalar particles, as shown in Fig. 6. In dimensional regularisation with d = 4 − 2ϵ, the box diagram involving a heavy vector of mass m V = √ x V m U , Fig. 6(a), where m U is the U 1 leptoquark mass, reads (omitting the coupling structure and gauge group indices) with u , v being the external massless spinors whose chirality is projected with P L = (1 − γ 5 )/2, and we defined m i = √ x i m U , with i = {1, 2, 3} as in Fig. 6. Similarly, in the case of the box diagram with a scalar exchange, Fig. 6(b,c), we distinguish two cases. The first case is when the scalar in question is the massless Higgs boson and the diagram reads where the same notation, as in the case of a box diagram with the vector exchange, applies. The second case corresponds to the exchange of the U 1 Goldstone modes, and since we work in the Feynman gauge, the diagram reads In the Feynman gauge, the couplings of the U 1 Goldstone modes, ϕ U , to the fermions are [26] where Q 1 L = q 3 L and Q 2 L = c q Q L − s q q 2 L , and similarly for leptons. In order to couple to the external Higgs bosons through couplings in Eq. (12), we need a chirality flip on the fermion lines corresponding to the masses m 1 and m 2 in Fig. 6(c), which reflects in the fact that the amplitude Eq. (A3) is proportional to these masses. On the other hand, the couplings of the G ′ (Z ′ ) Goldstone bosons, ϕ G ′ (ϕ Z ′ ) are with q ′ L = c q q 2 L + s q Q L , and similarly for leptons. Notice that they do not allow for couplings to q 3 L and ℓ 3 L . As we are interested in the Higgs-fermion operators with at least one third-generation fermion, we neglect the effects of the G ′ and Z ′ Goldstone bosons. In the following, we give explicit expressions for the loop functions that appear in the evaluation of Wilson coefficients, written in terms of f VL , f SL , and f SLR .

Triangle-diagram contributions to Higgs-quark operators
Another type of contribution to the Higgs-fermion operators comes from triangle diagrams as in Fig. 7. At zero momentum, these renormalize the top Yukawa coupling, while the next order in the momenta gives, after applying the equations of motion (EOMs) for the massless external fermions, a contribution to O  Before applying EOMs, these amplitudes match onto three independent structures (operators in the Green's basis). Following the notation from [58], these are The mapping of these onto the Warsaw basis operators reads where G i indicate the Wilson coefficients in Green's basis and we have used y t ∈ R. In order to compute the one-loop matching onto G uHD3 and G uHD4 , one needs to choose three different momentum configurations and solve the linear system where ⃗ G = (G uHD1 , G uHD3 , G uHD4 ) ⊺ , the matrix A EFT encodes the EFT amplitudes, and ⃗ T UV is determined from the triangle diagrams for different momentum configurations. The final result for the Wilson coefficients is momentum independent, and we can make a convenient choice for the three configurations (p 1 , p 2 , p 3 ) ∈ {(q, 0, −q), (q, −q, 0), (0, q, −q)}, such that the result can be written as where and , .

(A19)
In the following, we give explicit expressions for the loop functions that appear in the evaluation of Wilson coefficients, written in terms of f T and f T G where we define the function f T V G (x 1 , x 2 ; p 1 , p 2 , p 3 ) as

Finite contribution to Higgs-lepton operators from inverse see-saw mechanism
The Yukawa couplings of the U 1 Goldstone boson ϕ U in Eq. (A4) give contributions to the Higgs-lepton current operators O (1,3) Hℓ through the diagrams in Fig. 8.
FIG. 8. Diagrams contributing to the one-loop matching of the Higgs-lepton current operators when the inverse seesaw mechanism is implemented.
Expanding in 1/m R and keeping only the leading term gives the contribution to dimension-6 operators. Omitting Yukawa couplings, both diagrams give Notice that the Yukawa coupling involving S L in Eq. (A4) is proportional to m R , so the dependence on m R will cancel. Adding the two diagrams we get contributions to [O (1) Hℓ ] and [O Hℓ ].

Results: finite parts
In the following, we will present the results for the finite parts of the Wilson coefficients from the one-loop matching procedure. In the previous sections, we have shown the contributions in the full UV theory. In order to obtain the matching, both tree-level and one-loop contributions in the SMEFT must be taken into account. Diagrammatically, the UV amplitudes correspond to the ones in Fig. 6 and Fig. 7, while the SMEFT amplitudes to Fig. 9 and Fig. 10 for boxes and triangles, respectively. Equating the amplitudes in both theories, the finite pieces (FP) from matching read The relevant diagrams for the matching of C HD , both in the SMEFT and in the full 4321 model, are shown in Fig. 11. Similarly to the case of the triangle diagrams above, the four-point function receives contributions from four operators in the Green's basis [58]: Therefore, we need to choose four different momentum configurations for the external Higgs states, p 1 , p 2 , p 3 , and p 4 , set up the system for G H□ , G HD , G ′ HD , and G ′′ HD , and solve for G HD = C HD . We obtain In the limit of vanishing vector-like fermion masses, x Q,L → 0, the loop functions f i=3,4,5 go to 1, while f 2 (x V , 0) = f 1 (x V ), thus reproducing the results of Ref. [59,60].
where C is a vector with all dimension 6 Wilson coefficients, including flavor indices, and A is the anomalous dimension matrix, that can be read from [34][35][36]: The anomalous dimension matrix depends on µ through the running of the SM parameters. Notice that the impact of dimension 6 operators on the SM parameters running is a dimension 8 effect [61] so it can be neglected. The solution to this differential equation is where P is the µ-ordered product. If we neglect the running of the SM parameters, Eq. (B3) reads

Leading-log running in yt
The one-integral term of the expansion in Eq. (B3) corresponds to the leading-log running. In this work, the most important leading-log contribution to the operators relevant for the EW fit is the one induced by the top Yukawa (g s does not induce any leading-log running in the sector relevant for our model): where y 2 t is the average of y 2 t (µ), Substituting A t , we recover the leading-log formulas given in section III B 1. Running from 2.5 TeV to m t , and using DsixTools [61], we getȳ t ≈ 0.87, which is the value we use for y t in this work.
2. Next-to-leading-log running in yt and gs The two-integral term of (B3), or the log 2 term of (B4), give the next-to-leading-log running. Including y t and g s effects, using the tree-level expressions of the third-family Wilson coefficients of Table III, and running from the UV, we get: [C [C These expressions include the running of the SM parameters if for y t we use the averageȳ t of Eq. (B6), and for g s , the averageḡ Running from 2.5 TeV to m t , and using DsixTools [61], we getḡ s ≈ 1.1.

Leading-log running in gL
The weak running in g L gives contributions to C Hℓ and C Hq , but not to the corresponding singlet Wilson coefficients. Taking the leading log, neglecting the running of g L , and using that [C ℓℓ ] ijkl is real, we find that the third-family contributions are Hℓ ] 33 δ ij log Similarly, using that [C (1,3) qq ] ijkl is real, we find Hℓ ] 33 δ ij log where the sum on V = G ′ , Z ′ is understood to be over the G ′ and Z ′ contributions to C (1,3) qq . In terms of the model parameters, the Wilson coefficients read and similarly where The vertex modifications of the Z-boson are where δ U (T 3 , Q) is a family-universal contribution that is given by To be consistent with the one-loop matching, we include one-loop corrections in y t . In our model, they are given by the diagrams in Fig. 12. Notice that neglecting terms O(s q,ℓ ) and O(g 2 s,Y /g 2 4 ), diagrams (a) and (b) only affect the purely third-family vertices. In principle, diagram (c) affects the Z-mass, but due to the {α EM , m Z , G F } input scheme we are using this correction is translated into a correction of the flavor-universal shift δ U and the W -mass we show below. The one-loop corrections in y t to the Z-vertex modifications relevant for the EW fit (which excludes the Ztt coupling) in our model are If we neglect the y t running, these finite pieces exactly cancel with the non-log terms of the one-loop correction to Z →b L b L in Eq. (C17). Moreover, when we insert the four-fermion operators at the tree level (see Table III) in the one-loop corrections to all vertex modifications relevant for the EW fit, Z ′ and G ′ contributions systematically cancel, leaving only the leptoquark contributions. The variations of the W couplings relevant for the EW fit (which excludes the W tb vertex), including the one-loop corrections in y t (diagram (b) of Fig. 12) are ℓq ] 3333 δ i3 δ j3 , (C20) The W mass modification, including the one-loop corrections in our model, is (C22)

Electroweak observable variations
To construct the likelihood for the EW fit, we use the observables O of Tables 1 and 2 of [38], that collect the measurements presented in [62][63][64][65][66][67][68][69][70][71][72], and the SM predictions. Some sub-blocks of the observables are correlated. In Table VII and Table VIII we provide the correlations among these sub-blocks for the Z-pole observables [62] and the W -pole observables [67] respectively.  We calculate the variation of the observables, ∆O = O NP − O SM , at leading order in the EW boson vertex modifications and δm W , as defined in Eq. (C1), and neglecting the SM fermion masses. The expressions of the EW boson vertex modifications and δm W in terms of the SMEFT Wilson coefficients are given in Appendix C 1. To write the variation of the observables, it is convenient first to define where [O] SM represents the SM prediction of the observable O at the tree level. The variation of the Z-pole observables is then The variation of the W -pole observables, working in the limit V ij = δ ij , is It is not an observable of the EW fit itself, but can be calculated as a function of some of the EWPO [62,63]. Its variation, as function of the Z-vertex modifications is