Low-energy phenomenology of scalar leptoquarks at one-loop accuracy

We perform a complete study of the low-energy phenomenology of S1 and S3 leptoquarks, aimed at addressing the observed deviations in B-meson decays and the muon magnetic dipole moment. Leptoquark contributions to observables are computed at one-loop accuracy in an effective field theory approach, using the recently published complete one-loop matching of these leptoquarks to the Standard Model effective field theory. We present several scenarios, discussing in each case the preferred parameter space and the most relevant observables.


Introduction
The Standard Model (SM) of particle physics provides an excellent description of physical phenomena in a wide range of energies and scales. Despite no direct evidence for new physics emerged in direct searches at the LHC, for several years now some low energy measurements continue to show significant deviations from the respective SM predictions, which fuel the hope that some New Physics (NP) might be lurking somewhere at the TeV scale. The most significant and robust deviations, that we take into account in this work, are the following: • a deficiency in LFU ratios of rare B decays in muons vs. electrons, R(K ( * ) ) = Br(B → K ( * ) µµ)/Br(B → K ( * ) ee) [12][13][14][15], • deviations in differential angular distributions of the B → K * µ + µ − decay, as well as in several branching ratios of b → sµµ processes [16][17][18][19][20][21], • a longstanding deviation from the SM prediction in the muon anomalous magnetic moment (g − 2) µ [22,23].
While also other measurements show deviations from the SM, those above stand out and have been the focus of a large amount of theoretical and experimental effort. In all cases, large theory efforts for improving the SM predictions (often very challenging) have been undertaken, and several experimental endeavours and analyses have been set up for confirming, disproving, or providing cross-checks, for the anomalies. Indeed, new measurements scheduled to appear within the next few years are expected to clarify the nature of all these anomalies. A confirmation for the presence of new physics in any one of these observables would of course be revolutionary in our understanding of physics at the TeV scale. For the same reasons, an equally large effort has been put into finding possible new physics explanations. In case of the B anomalies, leptoquarks (LQ) at the TeV scale can provide good explanations, even combining neutral and charged-current anomalies. If they couple to both left and right-handed muons, also the muon anomalous magnetic moment could be addressed. In all scenarios, in order to find a good explanation it is necessary to consider the constraints imposed by a large set of observables generated both at tree-level and radiatively. In some cases, renormalization group evolution (RGE) of the operators generated at the matching scale down to the scale of the observables represent -2 -
In case of vector leptoquarks, such finite terms are calculable only in ultravioletcomplete models, thus making the analysis necessarily model-dependent (see e.g. [28][29][30][31] for analyses of specific gauge models of lepton-quark unification). On the other hand, scalar leptoquarks can be considered as self-consistent simplified models, and all observables can be computed precisely in terms of the LQ couplings and masses. A particularly promising set of LQ to address the observed anomalies are the S 1 = (3, 1, 1/3) and S 3 = (3, 3, 1/3) representations. 1 Several works have been dedicated to study their phenomenology. The S 1 leptoquark has been considered as possible mediator for all anomalies [32][33][34][35][36][37][38][39][40][41], with varying degree of success. S 3 , instead, has long been recognized to be a very good candidate to address the deviations in the b → sµµ transition [42][43][44][45][46][47][48][49]. Finally, the combination of both leptoquarks has been considered as a good combined explanation of charged and neutral-current B-anomalies [50,51] and possible ultraviolet (UV) completions have been proposed in terms of a composite Higgs model [52], combining flavour anomalies with a solution to the Higgs hierarchy problem, as well as in the framework of asymptotically safe quantum gravity [53].
More recently, one-loop computations of several observables in this model have been published [54][55][56][57]. The approach adopted in these works is to compute directly in the model the dominant one-loop contributions to the desired observables. This methodology is however prone to missing possible relevant effects, and is not suitable to be systematically generalizable.
In this work we aim to perform a complete one-loop analysis of the S 1 + S 3 model, focussed at addressing the anomalies listed above, while being consistent with all relevant experimental constraints. We adopt an approach based on effective field theories (EFT), leveraging on our previous work [58] where the complete one-loop matching of the S 1 + S 3 model to gauge-invariant dimension-six operators of the SMEFT, in the Warsaw basis [59], is presented. The EFT approach is designed to factorize the UV-dependent part of the problem, i.e. the UV matching, from the purely low-energy one. The latter involves RGE of the EFT coefficients to the energy scale of the observables and the computation of the observables at one-loop, within the EFT, see e.g. ref. [60] for a simpler case of a scalar singlet. As we shall describe, most of these steps are already available in the literature in complete generality. The complete one-loop UV matching, done manually as in [58], requires a substantial amount of work, however it is possible to proceed systematically without neglecting terms. Furthermore, this step is expected to become automatised in the near future. This will facilitate extending this work to include more observables, or to apply it to different UV models. In case of leptoquarks and low-energy observables, the use of EFT approaches is even more justified by the collider bounds from LHC, which put lower bounds on leptoquark masses close to the ≈ 1 TeV scale, see e.g. refs. [52,56] for recent reviews of pair production searches of S 1 and S 3 . Truncating the EFT expansion at dimension-six implies an implicit uncertainty in the evaluation of the NP contributions 1 We show the representations under the SM gauge group SU(3)c × SU(2)L × U(1)Y -3 -

JHEP01(2021)138
to observables, due to missing higher-dimension operators, that can be estimated being of O(E 2 /M 2 LQ ) or O(m 2 EW /M 2 LQ ) compared to the corresponding dimension-six contribution, where E is the typical energy of the process under consideration and m EW an electroweakscale mass. While former effects are completely negligible, the latter are ∼ 1% for TeV-scale leptoquarks, which do not affect the results in any sizeable way, given present day precision in the observables. 2 Our goal is to find interesting scenarios, within the S 1 +S 3 setup, capable of addressing one or more of the anomalies listed above, find the preferred region in parameter space, and discuss the most important experimental constraints in each case. Specifically, we first aim to quantify how well each leptoquark can address which set of anomalies, then we discuss combined explanations with both leptoquarks. Thanks to the complete one-loop matching, we also discuss limits on leptoquark couplings to the SM Higgs boson, arising from electroweak precision data and Higgs measurements.
In section 2 we present the S 1 + S 3 model, the methodology employed in the analysis, and present the list of all observables included in the fit, including a discussion of the relevant collider bounds, particularly those from Drell-Yan. The results for all scenarios considered are collected in section 3 and a discussion on future prospects can be found in section 4. We conclude in section 5. In appendix A we describe in details the LQ contributions to all the observables considered.

Setup
The Lagrangian for the two leptoquarks is the following where = iσ 2 , λ H1 , λ H3 , λ H3 ∈ R, (λ 1L ) iα , (λ 1R ) iα , (λ 3L ) iα , λ H13 ∈ C. We assume baryon and lepton number conservation 3 and we neglected quartic self-interactions between leptoquarks. The convention used for covariant derivatives is for a generic field Φ charged under the SM gauge group. We denote SM quark and lepton fields by q i , u i , d i , α , and e α , while the Higgs doublet is H. We adopt latin letters (i, j, k, . . . ) for quark flavor indices and greek letters (α, β, γ, . . . ) for lepton flavor indices. We work in the down-quark and charged-lepton mass eigenstate basis, where

3)
2 Dimension-eight terms could be relevant if they generate at tree-level an observable that is instead loop-induced at dimension-six. This, however, does not happen in this UV model for the observables under consideration. 3 See ref. [49] for an explicit setup forbidding baryon-violating couplings of S3 in a gauge model.

JHEP01(2021)138
and V is the CKM matrix. Except for the sign of gauge couplings, here and in the following we use the same notation specified in [58]. Integrating out at tree-level the two LQ, the following semileptonic operators are generated: The complete one-loop matching between the UV theory and the SMEFT in the Warsaw basis, as well as the definitions for the effective operators, are reported in [58].

Methodology
Our goal is to study the phenomenology of the S 1 + S 3 model described in the previous section, expressing the low-energy observables as functions of the UV parameters at oneloop level. Given the separation of scales between the LQ masses, assumed to be at the TeV scale, and the typical energy scales of the observables considered, the EFT approach is particularly suited for this goal. In fact, it allows to separate the complete procedure in a sequence of steps, which can be generalised to be applicable also to other UV scenarios. Going from the ultraviolet to the infrared, the matching procedure allows to pass physical thresholds, i.e. to integrate out heavy fields while defining a new EFT for that energy range, while the renormalization group evolution (RGE) allows to change the scale within an EFT approach. In our specific case, we have the following steps: JHEP01(2021)138 5. The expression of the low-energy observables and pseudo-observables in terms of the LEFT Wilson coefficients, taking into account contributions that arise at one-loop level within the LEFT, from the operators generated already at the tree-level. 4 By combining everything, we obtain expressions for the observables as a function of the parameters of the scalar leptoquark model at the TeV scale; in such a way, experimental bounds on low-energy data can be used to set constraints on the S 1,3 couplings. On the other hand, the intermediate steps provide model-independent expressions for observables in terms of EFT Wilson coefficients, which might be exploited in other NP scenarios. For a generic EFT coefficient we can separate a contribution arising at the tree-level from one arising at one-loop (2.5) Working at one-loop accuracy, the RGE, one-loop matching between SMEFT and LEFT, and the one-loop matrix elements to the observables, should only be considered for treelevel generated coefficients, C (0) i (in our case, those in eq. (2.4)). For the loop-generated coefficients, C (1) i , only the tree-level matching conditions from SMEFT to LEFT, and treelevel matrix elements should be included, the other contributions giving terms which are formally of two-loop order and that could be of the same order as neglected two-loop matching conditions.
The exception to this is in the RGE due to QCD from the TeV to the GeV scale, for example in four-quark operators contributing to ∆F = 2 observables. In this case the RGE contribution is well known to be important, also due to the large separation of scales, which gives to this effect a parametric enhancement with respect to the neglected two-loop corrections even if four-quark operators are generated at one-loop.

Observables
One of our main goals is to provide, with the S 1,3 model, a combined explanation for the hints of non LFU in the neutral and charged current semileptonic B-meson decays, namely to account for the experimental measurements of R K ( * ) and R D ( * ) , and of the deviation in the muon anomalous magnetic moment (g−2) µ . The leptoquark couplings involved in these observable enter also in the other low-energy observables (or pseudo observables), both at tree-level or one-loop level. Therefore, to quantify how the S 1,3 model can consistently explain the observed anomalies, one should take into account a set of low-energy data as complete as possible. In tables 1, 2, and 3, we show the list of low-energy observables that we analyze, together with their SM predictions and experimental bounds.
In appendix A, these low-energy observables are discussed in length. We will explicitly show, as functions of the parameters of the S 1,3 model, tree-level contributions together with dominant one-loop effects, while in the numerical analysis the full set of one-loop corrections is considered. Some of the considered observables vanish or are flavor-suppressed at tree-level, for example meson-mixing ∆F = 2 processes, τ → 3µ and τ → µγ LFV Br(K + → π + νν) 8.64 × 10 −11 [77] (11.0 ± 4.0) × 10 −11 [78] Br  interactions or τ → µφ(η, η ) decay; in such cases the inclusion of one-loop contributions is relevant and might bring non negligible changes in a global fit of the low-energy data. From the observables listed above, and their expression in terms of the parameters of the model, LQ couplings and masses, we build a global likelihood as: is the expression of the observable as function of the model parameters, µ i its experimental central value, and σ i the uncertainty. These are all discussed in appendix A. From the χ 2 built in this way, in each scenario considered we obtain the maximum likelihood point by minimizing the χ 2 , which we use to compute the ∆χ 2 ≡ χ 2 −χ 2 min . This allows us to obtain the 68, 95, and 99% CL regions. In the Standard Model limit we get a χ 2 SM = 101.0, for 50 observables.

JHEP01(2021)138
Observable Experimental bounds Z boson couplings Appendix A.12 scan over all the parameter space 5 and select only the points with a ∆χ 2 less than the one corresponding to 68 and 95%CL. The points obtained in this way also reproduce the CL regions in parameter space obtained by profiling. With this set of parameter-space points we can then plot any observable evaluated on them.

Collider constraints
Leptoquarks are also actively searched for at high-energy colliders. Their most important signatures can be classified in three categories: i) pair production, ii) resonant singleproduction, and iii) off-shell t-channel exchange in Drell-Yan processes, pp → + − or ν. See e.g. refs. [101,102] for reviews.
The pair production cross section is mostly independent on the LQ couplings to fermions, unless some are very large, and thus provides limits which depend only on the LQ mass and the branching ratios in the relevant search channels. We refer to refs. [37,52,56,102] for reviews of such searches. Once the branching ratios are taken into account, the most recent ATLAS and CMS searches using an integrated luminosity of ∼ 36fb −1 put a lower bound on the S 1 and S 3 masses at ≈ 1 TeV or less. At present, limits from single production are not competitive with those from pair production and Drell-Yan [102].
Leptoquarks can also be exchanged off-shell in the t-channel in Drell-Yan processes. The final states most relevant to our setup are ττ , τν, and µμ. The limits on LQ couplings as a function of their mass from neutral-current processes can be taken directly from [37,102,103] (see also [104][105][106][107] for other studies of dilepton tails in relation with B-anomalies) while the mono-tau channel in relation to the B-anomalies has been studied in [108][109][110][111][112][113] and at present it doesn't exclude the region of interest. Using the results from [37] we get the following 95% CL upper limits on the couplings relevant to our model for M 1,3 = 1 TeV, 5 For each numerical scan we collected O(10 4 ) benchmark points. For our more complex models (i.e. with up to ten parameters), this is quite demanding from the computational point of view; in order to efficiently scan the high-dimensional parameter spaces, we employ a Markov Chain Monte Carlo algorithm (Hastings-Metropolis) for the generation of trial points.
while the limits on λ 1L bτ and λ 1L bµ from the neutral-channel are in the non-perturbative regime since they induce only a |V cb | 2 -suppressed amplitude in cc → ττ , µμ. Since only the constraint on λ 1R cτ turns out to be relevant for our setup, we are justified taking the limits above derived one at a time.
We implement the constraints in eq. (2.7) in our global likelihoods by assuming that the signal cross section is dominated by the purely New Physics contribution, thus with a scaling ∝ |λ| 4 , and that it follows a Gaussian distribution around zero such that the 95%CL limits on single-couplings reproduce those in eq. (2.7), i.e. with standard deviation given by σ x = |λ lim x | 4 / χ 2 95%CL(1dof) , where for each coupling λ lim x is the limit in eq. (2.7) and χ 2 95%CL(1dof) ≈ 3.84: which is added to the global χ 2 from low-energy observables of eq. (2.6).

Scenarios and results
In this section we discuss several minimal models within the S 1 + S 3 setup, and how well (or bad) each of them is able to address the charged and/or neutral current anomalies, while remaining compatible with all the other experimental constraints. We denote the leptoquark couplings to fermions by: The main experimental anomalies driving the fit can be split in three categories: • CC : deviations in b → cτ ν transitions; • NC : deviations in b → sµµ transitions; • (g − 2) µ : deviation in the muon magnetic moment.
While our setup allows to keep all the above couplings in a completely general analysis, given the large number of parameters this would preclude a clear understanding of the physics underlying the fit. Furthermore, it can be interesting to consider only one leptoquark or to focus on one specific experimental anomaly. For these reasons we take a step-by-step approach by starting with single-leptoquark scenarios and switching on the couplings needed to fit a given set of anomalies. In all cases, we keep the complete likelihood -10 -

JHEP01(2021)138
Model The models are thus defined by the leptoquark content and the set of active couplings, which, for simplicity, we assume to be real. We have considered the models detailed in table 4, for each of which we allow the couplings listed in the third column of the table to be non-vanishing in our global fit. We first analyze single mediator models and study their potential to address as many anomalies as possible. In each case we point out the main tensions which prevent a combined explanation of all anomalies. Then, we move on to study models involving both leptoquarks. In the first we only allow left-handed couplings, λ 1L and λ 3L , as this possibility has better chances to find motivation in a scenario in which the flavor structure is determined by a flavor symmetry, see e.g. [51,52]. In the second we switch on also some of the S 1 couplings to right-handed fermions, and aim to provide a combined explanation for all three anomalies. Finally, we study the limits on the leptoquark potential couplings to the Higgs, which is an analysis largely independent on the couplings to fermions and requires to consider different observables than those studied in the main fit, see also [57]. In any given model there is, of course, no particular reason to expect the exact flavor structures implied by table 4. For instance, the couplings we set to zero will be radiatively generated. In our bottom-up approach we assume them to be small enough at the matching scale that the observables in the fit are not impacted in a sizeable way. In a more top-down approach one might have expectations on the size of these terms based on the UV picture, such as due to the presence of approximate flavor symmetries or other flavor-protection mechanisms.

JHEP01(2021)138
In the numerical analysis we fix for concreteness values of leptoquark masses equal to M 1 = M 3 = 1 TeV. While this is borderline with the exclusion limits from pair production, discussed previously, the results do not change qualitatively by increasing slightly the masses. Since most of the observables driving the fits scale as λ 2 /M 2 , with a good approximation this scaling can be used to adapt our fits to slightly larger masses. 6 We note that the future limits on LQ masses from HL-LHC are expected to not go much above 1.5 TeV [52].
Concerning our specific benchmarks, the choice of active couplings in each case is guided by some simple phenomenological observations (more details on each concrete model can be found in the relevant subsections below): 1. Since the observed deviations in B-decays involve LQ couplings to second and third generation, and given the strong constraints on s ↔ d quark flavor transitions, couplings to first generation of down quarks can only play a minor role in the fit of B-anomalies and are thus set to zero (note that even in this case, due to the CKM matrix, effects in up-quark observables are present, for instance D-meson mixing).
2. Hints to LFU violation in rare B-decays, combined with the deviations observed in B → K * µ + µ − , suggest that the LQ couplings to muons should be larger than those to electrons. We consider, for simplicity the case in which b → s anomalies are entirely explained by muon couplings and set to zero the couplings to electrons.
3. The S 1 couplings to µ R and c R or t R do not contribute to b → s , nor to b → c(τ / )ν, however are relevant for fitting the observed anomaly in the muon anomalous magnetic moment, which gets the main contribution from the couplings to b L µ L and t R µ R (cf. eq. (A.116)).
Details for all models are given in the following subsections.

Addressing CC anomalies
This LQ can address the deviations in R(D) and R(D * ) with only two couplings: λ 1L bτ and λ 1R cτ . They generate at tree-level a contribution to the semileptonic scalar and tensor operators C (1,3) lequ at the UV matching scale, eq. (2.4), which then run down to the GeV scale. The best fit region is entirely determined by the following few observables: The results from the fit, assuming real couplings and M 1 = 1 TeV, can be seen in figure 1. Since all the relevant low-energy observables scale with λ/M , the fit can be easily adapted to other masses. The left panel shows confidence level regions for the two couplings. The dashed lines are 95% CL constraints from single observables, and help illustrate the role of each observable within the global fit.

JHEP01(2021)138
Result from a fit in the two S 1 couplings λ 1L bτ and λ 1R cτ , for a leptoquark of 1 TeV. In the left panel we show the preferred regions in the plane of the two couplings, and the individual 2σ limits from the most relevant observables (except for (R(D), R(D * )), for which we show the 68%CL region). The black dot corresponds to the best-fit point. In the right panel we show where this preferred region is mapped in the plane of R(D) − R(D * ), together with the experimental combination from HFLAV [67].
In the right panel we show how the 68% and 95% CL region from the global fit of the left panel maps in the R(D ( * ) ) plane. This is almost degenerate, due to the approximate linear relationship between these two observables in the present model. We overlay as gray lines the CL ellipses (for 2 degrees of freedom) from the HFLAV global fit of the two observables [67] (specifically, the Spring 2019 update).
The model is successful in fitting the deviation in R(D ( * ) ) within the 68%CL level, with smaller values of R(D * ) (or larger R(D)) preferred by the fit. Improved measurements of R(D ( * ) ) can test this setup due to the precise linear relationship among the two modes predicted by the model, as well as improved Drell-Yan constraints.
The best-fit point, for M 1 = 1 TeV, is found for

Addressing NC anomalies
One may attempt to fit neutral-current anomalies in b → s from the one-loop contributions from S 1 . This scenario has been considered for the first time in [34]. Significant contributions to ∆C µ 9 may only come from the two muon couplings λ 1L bµ and λ 1L sµ , whereas the universal contribution is always negligible (C u 9 ≈ 0), cf. eqs. (A.6), (A.9). B s -mixing, cf. eq. (A.74), and B → K ( * ) νν, cf. eq. (A.49), put strong constraints on the product of the two couplings. Thanks to the different scaling of these observables and ∆C µ 9 on the leptoquark couplings, the limits can be avoided by a suitably large leptoquark mass, as can be seen in figure 2 (left). For M 1 3 TeV it is possible to avoid these limits while having couplings still in the perturbative range (see also [36,38]). Nevertheless, even while marginally evading the B s -mixing constraint, the ∆C µ 9 deviation remains in ≈ 2σ TeV. In both plots the green region is the 1σ favourite one from ∆C µ 9 .
tension with the bound on λ 1L bµ arising from the LFU limit from τ decays, |g µ /g e | (see appendix A.8). The situation is illustrated in the right panel of figure 2. This is slightly exacerbated by a ∼ 1σ deviation in the opposite direction measured in |g µ /g e |.
We thus conclude that the S 1 leptoquark is not able to fit neutral-current anomalies while remaining completely consistent with all other constraints. The situation regarding NC anomalies is not modified significantly by letting also the other couplings vary in the fit. This issue could be avoided by allowing a mild cancellation by tuning a further contribution to this observable, possibly arising from some other state. Fixing M 1 = 6 TeV we find the best-fit point for:

Addressing (g − 2) µ,e
Thanks to the m t enhancement of the left-right contribution, cf. eq. (A.111), S 1 is a good candidate to address the observed deviation in the muon anomalous magnetic moment. The leading contribution is given numerically by (see eq. (A.116)): The observed deviation can thus be addressed for small couplings, and no other observable is influenced significantly. Analogously, it is possible to address the (smaller and less significant) deviation in the electron magnetic moment, see table 2.
A combined explanations of both deviations with a single mediator was thought not to be viable, due to the very strong constraint from µ → eγ, see e.g. refs. [114]. More recently, in the updated versions of [115,116] it was realized that a possible way out is to align the S 1 leptoquark-muon couplings to the top quark, while the leptoquark-electron ones to the -14 -JHEP01(2021)138 charm. 7 In our formalism this can be achieved aligning the S 1 couplings as λ 1L : On the one hand, the strong limit from µ → eγ implies that the alignment described above must be held with high accuracy. On the other hand, radiative corrections to leptoquark couplings from SM Yukawas are expected to induce deviations from it. For this reason we will not investigate this direction further.

Addressing CC and (g − 2) µ
Following the previous subsections, the next natural step is to attempt to fit both R(D ( * ) ) and a µ anomalies with tµ . This setup is not simply the combination of those discussed previously. Indeed, due to the λ 1L bτ and λ 1R cτ couplings on the one hand (required to fit the charged-current B-anomaly), and λ 1R tµ , λ 1L bµ on the other hand (necessary to fit a µ ), sizeable contributions to τ → µγ are generated at one-loop, see eq. (A. 113). The values of λ 1L bτ,bµ and λ 1R cτ,tµ required to fit R(D ( * ) ) and (g − 2) µ would induce a too large contribution to this LFV decay. However, we find that the large contribution to τ → µγ arising from the product of λ 1L bµ λ 1R cτ can be mostly cancelled by the λ 1L bµ λ 1R tτ term, for λ 1R tτ ∼ 0.003, corresponding to a small tuning of approximately one part in 5. Such a small value does not affect any other observable in our fit. It should be noted that the only effect of this coupling is to tune this observable.
We show the preferred region in parameters space, obtained from our fit, in figure 3. In the upper two panels we show the fit in the (λ 1L bµ , λ 1R tµ ) and (λ 1L bµ , λ 1R tτ ) planes, including the preferred region from the global fit as well as the individual 95%CL limits from single observables, to help illustrating the physics behind the analysis. It can be noted that the observed value of the muon magnetic moment can be reproduced, and that τ → µγ requires a non-zero value of λ 1R tτ , as discussed above. Regarding the top-right panel, the single-observable constraint from τ → µγ (red dashed line) is shown by imposing that a µ is fixed to the value we get at the best-fit point. In the lower-left panel we show the fit in the couplings contributing to R(D ( * ) ), (λ 1L bτ , λ 1R cτ ), and how this preferred region is mapped in the plane of R(D) − R(D * ) (lower-right panel). Comparing to the allowed region in the same plane in the model studied in figure 1, we see that τ → µγ strongly reduced the allowed region, while still not preventing a good fit of the charged-current B-anomalies. Due to this reduction of the allowed (λ 1L bτ , λ 1R cτ ) parameter space, specifically with smaller values of λ 1L bτ , the points in the R(D * ) − R(D) plane line up more closely in a line than what is observed in figure 1. 7 We thank the authors of [115] and [116] for pointing this out to us. We conclude that the R(D ( * ) ) and (g − 2) µ anomalies can be addressed by the S 1 leptoquark, for perturbative couplings and TeV-scale mass. We find the best-fit point, (3.6)

Single-leptoquark S 3
We move on to examine S 3 , and we attempt directly a combined explanation of charged and neutral current anomalies. It is well known that S 3 provides a simple and good explanation for the deviations observed in b → s , thanks to its tree-level contribution to the partonic process. The couplings required are λ 3L bµ -λ 3L sµ , with small enough values that other observables do not pose relevant constraints. The leading contribution to both R(D ( * ) ) and C u 9 , instead, arises via the λ 3L bτ -λ 3L sτ couplings. For concreteness we fix M 3 = 1 TeV, but the fit would be very similar for a slightly larger mass.
Our results can be seen in figure 4. As expected, the model is successful in fitting ∆C µ 9 . The couplings to the tau allow to also fit C u 9 , while charged-current anomalies cannot be reproduced. The main limiting observables are B s -mixing and B → K ( * ) νν, as can be seen from the top-right panel.
The best-fit point, for M 3 = 1 TeV, is found for

S 1 + S 3 with LH couplings only
Models involving S 1 and S 3 with left-handed couplings have been first considered in [50][51][52]. In particular, in [51,52] it was shown how this setup could fit both charged-and -17 - neutral-current anomalies with couplings compatible with a minimally broken U(2) 5 flavor symmetry, albeit with a tension between R(D ( * ) ) and the B s -mixing constraint. Since then, new experimental updates on R(D ( * ) ) pushed the preferred region closer to the SM, thus also alleviating the tension with meson mixing. Here, we update the fit for this scenario, without assuming a priori a specific flavor structure for the relevant couplings.
The relevant couplings are λ 1L [bs]τ , λ 3L [bs]µ , and λ 3L [bs]τ . A first qualitative understanding of the model can be obtained by noticing the main roles of the various couplings with regard to the anomalies: In this model, the relative deviation in R(D) and R(D * ) from the respective SM values is predicted to be exactly the same, since it is only due to the same left-handed vector-vector operator generated in the SM. The contribution is given by the combination in eq. (A.29).
The most salient features of the fit are summarized in figure 5. In the top two panels we show the preferred regions in the λ 1L bτ − λ 1L sτ and λ 3L bτ − λ 3L sτ planes, together with the single-observable 2σ limits obtained fixing the other couplings to the global best-fit value. The favoured region in the λ 3L bµ − λ 3L sµ plane is very similar to the one of model S (figure 4 top-left), thus we do not show it again. The constraint from B → K ( * ) νν is avoided thanks to a slight cancellation between the tree-level contributions of the two leptoquarks [51], see eq. (A.43). There is a (small) leftover tension in the R(D ( * ) ) fit, due to constraints from B s -mixing. It should be noted that this tension grows with larger LQ masses (thus larger required couplings) since the deviation in R(D ( * ) ) scales as λ 2 /M 2 while the contribution to meson mixing goes as λ 4 /M 2 .
We also point out that the parameter-region preferred by the fit is compatible with the relations between couplings predicted by a minimally-broken U(2) 5 flavor symmetry, [51] and references therein. The case with |c U(2) | = 1 is shown with grey dashed lines in the upper panels.
In the lower two panels we show how the preferred regions in parameter space maps into the anomalous B-decay observables. As can be seen, it is possible to reproduce them within 1σ. The best-fit point, for M 1 = M 3 = 1 TeV, is found for (3.9)

S 1 + S 3 addressing CC, NC, and (g − 2) µ
From the previous sections it is clear that in order to address all anomalies, both S 1 and S 3 leptoquarks are required. NC anomalies are addressed only by S 3 , the muon anomalous magnetic moment only by S 1 , while R(D ( * ) ) receives sizeable contributions from both. For our most general analysis we keep ten active couplings: The results of our fit are shown in figure 6.
In the first row of figure 6 we show the preferred regions for the couplings relevant for the a µ fit. The situation is very similar to what already discussed for model S The couplings relevant for the R(D ( * ) ) fit are shown in the second row. They show a behavior very similar to the one already seen in the models S (CC+aµ) 1 The main contribution is due to the scalar+tensor operators generated via the λ 1R cτ λ 1L bτ couplings, but a sizeable contribution, which helps to improve the fit with respect to model S . Contrary to that case, however, here the preferred region avoids any tension with both B s -mixing and B → K ( * ) νν.
We do not show in figure 6 the preferred values for λ 3L [s,b]µ , which are necessary to fit ∆C µ 9 , since they are analogous to what we saw for model S (see figure 4 top-left). -20 -

Leptoquark potential couplings
In this section we study available constraints for the potential couplings of leptoquark with the Higgs boson in the third and fourth lines of eq. (2.1). There are four such couplings: λ H1 , λ H3 , λ H3 , and λ H13 . All contribute only at one-loop level in the matching to SMEFT operators, therefore possible phenomenological effects are suppressed both by a loop factor and by the LQ mass scale. We focus on effects of these couplings which are independent on the LQ couplings to fermions. We thus need precisely measured quantities in the bosonic sector of the SM. Obvious candidates are the gauge-boson oblique corrections measured at LEP [121]:Ŝ, T , Y , W , as well as the analogous effect for QCD, Z. All these parameters are measured at the per-mille level, and are able to constrain multi-TeV scale physics. Given the expressions in the Warsaw basis of [122] and our one-loop matching of the SMEFT to the LQ model, we obtain (see appendix A.13 for details) where in the numerical expressions for simplicity we fixed M 1 = M 3 = m TeV. The contributions to Y , W , and Z are instead at, or below, the 10 −6 level and thus completely negligible given the present experimental precision. The constraints on S and T from [117] are reported in  been also studied in [40], albeit not in the EFT approach. We checked that we agree once the EFT limit is taken into account.
The LQ couplings to the Higgs also generate at one-loop contributions to hgg, hγγ, and hZγ couplings. Since these are also loop-generated in the SM, the percent-level precision presently available for the Higgs couplings to photons and gluons couplings allows to probe heavy new physics. Loop contributions to other couplings, which arise at tree-level in the SM, are instead too small to have a sizeable impact. We thus consider the combined fit of Higgs couplings in the κ-framework where only κ γ and κ g are left free, and a constraint on σ/σ SM (Zγ) = κ 2 g κ 2 Zγ , which is however still not precisely measured, see table 5. The approximate contributions to these parameters in our model are given by (details in appendix A.13) Analogously to what presented above for flavour observables, we combine Higgs couplings and oblique constraints in a global likelihood. From this we find the maximum likelihood point and construct the 68, 95, and 99% CL regions in planes of two couplings, where the other two are marginalised. The results in the (λ H1 ,λ H3 ) and (λ H13 ,λ H3 ) planes are shown in figure 7 for M 1 = M 3 = 1 TeV. We observe that a limit of about 1.5 can be put on both λ H13 and λ H3 (right panel). This comes mainly from the contribution to theT parameter, eq. (3.12), which is quadratic in the two couplings and thus allows to constrain both at the same time. The λ H1 and λ H3 couplings, instead, are constrained mainly from their contribution to the hγγ and hgg couplings, eq. (3.13). We see that with JHEP01(2021)138 present experimental accuracy the limits are still rather weak, and there is an approximate flat direction which doesn't allow to put any relevant bound on λ H1 . This situation will marginally improve with the more precise Higgs measurements from HL-LHC [123]. The future expected 95%CL contours are shown as dashed blue lines. This however has no appreciable effect on the limits shown in the right panel, since those are dominated by the constraint on the T parameter, which will instead improve substantially from measurements on the Z pole at FCC-ee. A more detailed analysis of FCC prospects are however beyond the scope of this paper.

Comparing with literature
In recent months the S 1 + S 3 model at one-loop accuracy has been studied in refs. [55,56] for what regards the flavour anomalies, while ref. [57] studied electroweak and Higgs limits on the leptoquark-Higgs couplings. Given the similarity of the goals with out work, we discuss in this section the main differences. The most important lies in the approach used to calculate radiative leptoquark contributions to observables. While previous works employed direct computations of leptoquark loop contributions to the desired low-energy amplitudes, in this work we use an EFT approach, whereby the only model-dependent part of the computation is the one-loop matching to the SMEFT. As argued in the introduction, we believe such an approach has several advantages, the most important being the automatic inclusion of all new physics effects to all observables at leading order in 1/M 2 LQ expansion and to one-loop accuracy: there is indeed no need to simplify the computation neglecting given terms or couplings, for example all electroweak corrections are included automatically in our computation. 8 This approach thus allow us to study a larger set of observables than the ones considered in previous studies. One example is D 0 −D 0 mixing, which was not considered in [55,56]. Considering the benchmark points selected by both works, we find that all of them are excluded (by a large amount) by D-meson mixing, due to large λ 1(3)L 23

couplings.
In the benchmarks of [56] we also find a large tension in τ → µγ, while we are in agreement with the analytic expressions for the decay. We also find some tension in the leptonic decay D s → τ ν for several benchmarks points of both studies.
We should also point out that the Higgs-S 3 interaction proportional to λ H3 , cf. eq. (2.1), was not included in the analysis of ref. [57].

Prospects
In this section, we discuss the implications of future Belle II measurements of • LFV B decays induced at parton level by b → sτ µ (see appendix A.6); • B decays induced at parton level by b → sτ τ (see appendix A.6).
These processes, in fact, are particularly interesting for leptoquark scenarios aiming at addressing both neutral and charged-current B-anomalies. Both are induced at tree-level JHEP01(2021)138 Table 6. Future Belle II and LHCb sensitivities, at 95% C.L., for b → sτ µ and b → sτ τ observables. by S 3 and, by SU(2) L relations, the b → cτ ν τ transition, tree-level in the SM, is related to the FCNC transition b → sτ τ . Also, LFV is a natural consequence of leptoquark couplings once also the coupling to muons is considered, as required by neutral-current anomalies. While the LFV B-meson decays are already included in the global fits described in the previous sections, the current bounds on b → sτ τ observables are of the order of ∼ 10 −3 and thus too weak to set constraints on the model parameters. However, Belle II, with 50ab −1 of luminosity, will strongly improve the sensitivity, in particular for the branching fraction of the semileptonic decays. On the other hand, the Upgrade II of LHCb will set competitive bounds on the leptonic decay B s → τ τ . The relevant future expected limits at 95% C.L. for Belle II [124] and LHCb [125] are summarised in table 6.
In figure 8 we show how the preferred parameter-space regions for the models S 1 +S decay B + → K + τ µ and the decay B + → K + τ τ (normalised to the SM value). 9 The red horizontal lines correspond to the Belle II future bounds at 95% C.L. on Br(B + → K + τ + τ − ) at 5ab −1 (dashed lines) and 50ab −1 (solid lines), while the vertical ones represent the Belle II 50ab −1 prospect for Br(B + → K + τ µ). One can see that, in both scenarios, the predictions for both the non-LFV and LFV semileptonic B decay into τ are in the ballpark of the future Belle II sensitivity at 50ab −1 , while the expected bounds at 5ab −1 are still too weak to set significant constraints on the models. Furthermore, one can notice that the future measurements of b → sτ τ observables are constraining more strongly the parameter space of the S 1 + S 3 (LH) model than the one of the S 1

Conclusions
In this work we examined in detail and at one-loop accuracy the phenomenology of Standard Model extensions involving the two leptoquarks S 1 and S 3 , motivated by the experimental discrepancies observed in B-meson decays and in the muon anomalous magnetic moment (g − 2) µ .
To this aim, we performed global fits for several benchmark models to a comprehensive list of flavor and electroweak precision observables, each computed at one-loop accuracy, leveraging on our previous work [58]. For each model, we identify best-fit regions and major sources of tension, when present, and also provide prospects for B-decays to τ τ and τ µ, in the experimental scope of Belle II and LHCb.
It is found that models involving only the S 1 leptoquark can consistently address R(D ( * ) ) and (g − 2) µ anomalies, while a fully-satisfactory solution for b → sµµ anomalies is prevented by the combination of constraints from B s -mixing and LFU in τ decays. Conversely, the S 3 leptoquark when taken alone can only address neutral-current B-meson anomalies. A model with both S 1 and S 3 , and only left-handed couplings for S 1 , can address both B-anomalies but not the muon magnetic moment. Finally, allowing for right handed S 1 couplings makes it possible to fit also (g − 2) µ . Concerning the prospects for both the LF conserving branching fraction Br(B → Kτ τ ) and the LFV one Br(B → Kτ µ), they are found to be in the ballpark of the future expected sensitivity of Belle-II and LHCb.
A quick glance summary of the various models is provided by figure 9 where we show, for the best-fit point of each model, the deviations from the SM prediction of each of the most relevant observables studied in the global fit. The black dots and intervals represent 9 It should be noted that at tree-level in our model this ratio is the same for all decays involving the b → sτ τ transition, e.g. Bs → τ τ (see appendix A.6). Figure 9. For the best-fit points in each model studied in this work, we show the relative deviations from the Standard Model prediction in all observables, in terms of number of sigmas given by the experimental precision in that observable. The black intervals represent the experimental measurements the experimentally preferred values and uncertainties, and for each observable we normalize the x-axis to the corresponding uncertainty (i.e. we count the number of standard deviations). Detailed informations for each model can be found in subsections 3.1-3.4. A separate analysis is provided in section 3.5 for Higgs physics observables and electroweak oblique corrections, which put constraints on leptoquark-Higgs couplings.
To conclude, we find that the combination of S 1 and S 3 provides a good combined explanation of several experimental anomalies: charged and neutral-current B-meson anomalies as well as the muon magnetic moment. Their mass is necessarily close to the 1 TeV scale, particularly to address charged-current anomalies R(D ( * ) ), and is thus in the region that could still show some signals at HL-LHC, if they are light enough, but that will definitely be tested at future hadron colliders.

JHEP01(2021)138
In the next few years, several experiments are expected to provide concluding answers as to the nature of all these puzzles. While at this time it is still very possible that some, or all, of these will turn out to be only statistical fluctuations and will be shown to be compatible with SM predictions, the possibility that even only one will instead be confirmed is real. Such an event would have profound and revolutionary implications for our understanding of Nature at the smallest scales. The scalar leptoquarks considered here are very good candidates for combined explanations and could thus be the heralds of a new physics sector lying at the TeV scale.
• Any other notation/convention agrees with those employed in our work [58].
Finally, we give, for convenience, our numerical expressions in terms of:

A.1 Observables for b → s
One of the goals of the present analysis is to provide an explanation for the hints of non LFU in the neutral-current semileptonic decay of B-meson into K ( * ) . We study here the observables related to b → s + − , focusing on those for which it is possible to obtain a robust SM prediction, i.e. the LFU ratios R K and R K * in several q 2 bins, as well as the leptonic decay B s → µ + µ − . Their SM predictions and experimental measurements are reported in the first five lines of table 1. The standard notation for the effective operators relevant to these processes is where the C i are evaluated at m b scale and the operators O i are defined as The expressions for the Wilson coefficients of the above operators in terms of the LEFT ones are of the parameter space in which the tree-level term turns out to be small. The leading contributions are the ones proportional to the squared top Yukawa and to the large logs associated to the electromagnetic RG in the LEFT, but there might be relevant effects also from terms in which the product of four LQ couplings enters. Therefore, an approximated expression for these semileptonic Wilson coefficients, at m b scale, in terms of the model parameters is where η ij kl ≡ V * ki V lj and for S 3 we kept only the tree-level contribution and the electromagnetic universal and vector-like loop correction.
Therefore, the coefficients in eq. (A.5) are all vanishing at first approximation, apart from C sb 9(10) which can be obtained, using eq. (A.6), as The relevant observables for these operators are the LFU ratios R K [12,14,129] and R K * [13,15], the branching ratio of B s → µ + µ − [19,20,130], the angular observable P 5 in B → K * µ + µ − [21], and branching ratios in other b → sµµ transitions [18,[131][132][133]. In our analysis, we use the results from the global fit of b → s observables done in [66]. In particular, with negligible S 1 couplings to right-handed muon λ 1R iµ ≈ 0, the fit relevant for us is the one with a LFU-violating contribution to muons along the left-handed combination ∆C sbµµ 9 ≡ C sbµµ 9 − C univ 9 ≈ −C sbµµ 10 , and a flavor-universal contribution along the vector-like direction C univ with a correlation of ρ ≈ −0.5. 10 The numerical expressions we get in our setup are, assuming real LQ couplings: 10 We use the results updated in April 2020 from P. Stangl's slides at https://conference.ippp.dur.ac.uk/ event/876/timetable/.

JHEP01(2021)138
We report for completeness the LFU ratios R K and R K * dependence on these Wilson coefficients: where the SM prediction is known up to O(1%) electromagnetic corrections [134] and we neglected subleading contributions due to quadratic terms or imaginary parts of the Wilson coefficients

A.2 Observables for b → cτ ν
In this section, we analyze the effects of the S 1,3 model on the observables related to b → c ν, in order to account for the hints of non LFU in the charged-current semileptonic B decays. In the following, we will at first assume that the NP contribution to the τ ν τ channel is the dominant one and all the others might be neglected. Indeed, as shown in the next section, LFU in light leptons has been checked at the percent level, this limits the size of the possible effect in the muon channel in R(D ( * ) ) to be sub-leading.
The relevant four-fermion effective Hamiltonian, at the µ = m b scale is where the C i are evaluated at m b scale and the operators O i are defined as They coincide with the ν τ τ bc matrix elements of the hermitian conjugates of the LEFT O νedu operators. The expressions for the Wilson coefficients in terms of the LEFT ones are given by approach we also include finite corrections arising from the one-loop matching between the LEFT and SMEFT [64], as well as in the SMEFT to LQ matching. These effects have a comparable impact to the RGE, as we will show below. Going from LEFT coefficients at the m b scale up to SMEFT ones at 1 TeV we get: where in the last expression we neglected small CKM and loop-suppressed terms propor- , which are however kept in the numerical analysis. The tree-level matching between SMEFT and the LQ theory is given in eq. (2.4), which sets the relation [C In terms of low-energy coefficients at the m b scale, we obtain the relation C S L ≈ −7.43C T . The observables mostly relevant to the b → cτ ν transition, for which a robust SM prediction is possible, are The approximate dependence on the non-vanishing EFT coefficients in eq. (A.18), valid up -31 -

A.3 LFU in b → cµ(e)ν
Relaxing the hypothesis that NP is coupled only to the third generation, we can have also corrections to b → cνe and b → cνµ interactions. Lepton Flavour Universality between -32 -

JHEP01(2021)138
muons and electrons in these charged-current transitions has been checked experimentally at the percent level [71,72]. Since these processes are also used to extract the CKM element V cb , if also new physics affects them a careful analysis should be performed in general, see e.g. [137]. In our case, since we assume negligible leptoquark couplings to electrons, we avoid this issue by using only the LFU ratio to constrain new physics. Since this ratio is insensitive to V cb , and a global fit for this SM parameter gets contributions from a large number of observables, we argue that the extraction of the CKM element wouldn't change in a sizeable way by removing this single observable. This ratio has been measured by both Babar [71] and Belle [72] collaborations: where the first uncertainty is statistical and the second systematic and in deriving this for Babar we took into account the correlation between systematic uncertainty in the muon and electron channels. Combining these two results, we get The decay rate into light leptons, as function of the same operators as in eq. (A.13) (but with = µ, e instead of the tau) is given by [138] Br where C X L+R ≡ C X L + C X R . Assuming negligible couplings to electrons, the leading contributions to the LFU ratio are given by where the expression for 2C µ V L is the same as in eq. (A.29) with τ → µ.

A.4 D s → τ ν
In the LQ model the leading contributions to this leptonic decay arise at the tree-level. Given the limited present experimental accuracy, and the fact that the tree-level contribution is expected to be the dominant one in these modes, we limit ourselves to this. The branching ratio of the D s → τ ν decay can be expressed as a function of the LEFT coefficients L V,LL νedu and L S,RR νedu , evaluated at the charm quark mass scale, The tree-level matching to the LQ model, augmented by the QCD RG evolution, gives The dominant contribution arises from the interference with the SM, thus from the τ channel: where Br(D + s → τ + ν τ ) SM = (5.169 ± 0.004) × 10 −2 [73]. The experimental measurement for the D + s → τ + ν τ branching ratio is [74] Br(D s → τ ν) A.5 B → K ( * ) νν, K + → π + νν, and K L → π 0 νν The relevant parton level processes d i → d j ν α ν β are described by the effective four-fermion Lagrangian with the following ∆F = 1 operators The relations between the Wilson coefficients of the operators above and the LEFT ones are given by In our model, at tree-level only C L is not vanishing. While at loop level a contribution to C R is generated, it is suppressed by small down-quarks Yukawa couplings, and thus it is completely negligible. The leading contributions to the Wilson coefficients at m b scale, in terms of the UV parameters, are where M ∼ M 1 , M 3 is the matching scale to the UV model. For B → K ( * ) νν it is customary to define the ratio with the (clean) SM prediction, The experimental limits at 95% C.L. are [76] These ratios can be expressed in terms of 2 real parameters > 0 and η ∈ [−1/2, 1/2] [75]: (A.47) The κ η parameter depends on form factors and its numerical value in [0,q 2 max ] is 1.34 ± 0.04, where q max is the kinematic limit of 22.9 GeV for K and 19.2 GeV for K * . At the leading order, in our model, C R = 0 and therefore η αβ = 0 and R ν K = R ν K * = 2 . The SM Wilson coefficient is given by The bounds of eq. (A.45) can probe a region of the parameter space in which the B → K ( * ) νν cross section is dominated by the squared BSM amplitude. Assuming real couplings we can provide approximate numerical expressions: where m i = M i /TeV.

A.6 b → d i τ τ and b → d i τ µ decays
We consider the leptonic branching ratios: 56) and the semileptonic ones: At partonic level, these processes are induced by the same low-energy operators as b → s , see eqs. (A.5), (A.6), with obvious substitutions for the flavor indices. We limit ourselves to C 9,10 contributions; furthermore we restrict to a tree-level analysis for the τ τ modes, whose current experimental bounds are too weak to be included in our global fit. In our model, the tree-level generated effective vertices lead to left-handed four-fermion interactions with C 9 sbτ τ = −C 10 sbτ τ , which are the ones that interfere with the SM processes. Concerning the τ modes, we have: where C 9,10 are defined as in eqs. (A.5), (A.6), C (9−10) ≡ C 9 − C 10 and we have Numerically, we get: (A.63) The present 95% CL limit on this branching ratio is 6.8 × 10 −3 [127], so almost a factor ∼ 10 4 larger than the SM prediction. The lepton flavor violating modes are zero in the SM and we need the full branching ratio expressions. Neglecting the muon mass (amounting to an O(10%) effect), the B d i → τ − µ + branching ratios can be expressed as [141] Br where the input parameters are f Bs = 0.224GeV [73], m Bs = 5.37GeV and τ Bs = 1.51 × 10 −12 s, f B d = 0.190GeV m B d = 5.279GeV and τ B d = 1.5 × 10 −12 s [74]. The Br(B d i → µ − τ + ) case is analogous, via exchange of lepton indices. In our analysis we keep m µ effects as well as the complete one-loop matching contributions. Summing the τ − µ + and µ − τ + modes, we get these approximate numerical expressions: Finally, for the LFV three body decays Br(B + → K + τ ± µ ∓ ), we get the integrated expressions for the branching ratios for B → K ( * ) τ µ from [141]. The leading contribution is given by where the omitted terms that do not contribute to ∆F = 2 processes and the ones proportional to lepton Yukawa couplings. We also defined In our fit we use the limits on ∆F = 2 operators from the UTFit collaboration [83]. 12 The bounds on the ∆F = 2 Wilson coefficients (or their real and imaginary parts, where appropriate) are given in terms of a scale Λ NMFV , defined via C = F The coefficients are evaluated at the same scale Λ NMFV , which being close to the TeV in most cases is close enough to the LQ masses we are interested in that we can neglect the residual running to the LQ matching scale. All constraints are collected in table 2. In terms of SMEFT Wilson coefficients, and LQ couplings, the relevant non-zero combinations are 12 We employ the results updated at the 2018 La Thuile conference https://agenda.infn.it/event/ 14377/timetable/.
-39 -JHEP01(2021)138 We remark that the bounds in ref. [84] assume a single non-vanishing coefficient at a time, so that the latter bounds could in principle be evaded by fine-tuning C 1 D and C 4 D ; we neglect this possibility, also because in the scenarios we focus on, the leading constraint is the one from B s -mixing.

A.8 LFU in τ decays
Lepton flavor universality has been tested in τ decays at the permille level, in particular the strongest constraints are obtained from the ratios [85] The decay amplitudes can be described by four-lepton operators at the tau mass scale, which are one-loop suppressed in our model. Given this suppression and the experimental precision, we can confine ourselves to contributions interfering with the SM one, which comes from a (V − A) 2 chiral structure; we can thus write: In terms of these ratios, the LFU quantities in eq. (A.75) are given by: In the LEFT these operators do not run at the leading order in the EFT expansion. The matching to the SMEFT at a scale µ W is given by The leading terms in the one-loop matching from the LQ model to the SMEFT coefficients [58] give: From the complete expression of the matching, and neglecting at the end LQ couplings to electrons and to second generation left-handed quarks (which are expected to be suppressed due to quark flavor constraints), the leading numerical contributions are: (A.82)

A.9 LFV decay τ → µφ
The relevant effective operators for the τ → µss transition in the LEFT are: At the tree-level, in our model the only non-vanishing contribution is Since the coupling to the s quark might be flavor-suppressed, loop contributions from LQ couplings to third generation quarks can potentially give a sizeable contribution. We evaluate the Wilson coefficients at the τ lepton mass scale, including all contributions up to one-loop arising from the matching to the UV model, the RGE, and the matching between SMEFT and LEFT. Simplified analytical formulae, which can help to understand the physics underlying this observable, can be derived by putting to zero electroweak corrections (except for the -41 -

JHEP01(2021)138
In principle we should also include one-loop contributions in the LEFT, from operators generated by integrating out the leptoquarks at tree-level, to the expression above in order to consistently obtain finite corrections. This might require calculation of non-local hadronic form factors, which is beyond the purpose of this work. However, one can realize that these are necessarily finite QED corrections, thus very suppressed with respect to the logarithmic correction from the QED RGE, which has a factor ∼ log M 2 LQ /m 2 τ ∼ 13 for a 1 TeV leptoquark. For this reason we neglect them. Taking all couplings real we can obtain the approximate expression: where we normalize the LH couplings to the strange quark with |V ts | ≈ 0.041 and the last line is dominated by the tree-level contribution. The experimental limit on this branching ratio at 95% C.L. is [74] Br(τ → µφ) < 1.00 × 10 −7 . (A.93)

A.10 LFV decay τ → 3µ
New physics contributions to the three body LFV decay τ → 3µ can be described by the effective Lagrangian at the tau energy scale: The branching ratio is given by [143,144] Br(τ → 3µ) = m 5 In terms of the LEFT operator basis in appendix B.2, the EFT coefficients above are given by: where in the second equalities we took into account the symmetry of flavor indices in the LEFT coefficients. Note also that the overall sign in X γ depends on the convention used for covariant derivatives, our is the same as [143]. In our LQ model, all these operators arise at loop level, except the scalar ones Q S,LL(RR) which are not generated at dimension-six in the SMEFT and therefore we will not consider in the following. There are three different types of contributions: the SMEFT-LEFT matching at tree-level, involving SMEFT operators that in the model are generated at one-loop level; the one-loop SMEFT-LEFT matching induced by tree-level four-fermion semileptonic operators; one-loop contributions in the LEFT, generated by penguin diagrams from the tree-level generated operators O V,LL ed and O eu . The latter are purely QED loops, so the only phenomenologically relevant contribution is the log-enhanced one, corresponding to the QED-induced RG evolution. Taking into account these different contributions and retaining the terms that should be dominant in most of the relevant region of parameter space, which is to say the ∼ λ 4 , ∼ y 2 t λ 2 ones and the LEFT RG running involving third quark generation, we can express the LEFT Wilson coefficients, evaluated at τ mass scale, as a function of the UV model parameters: In most cases, unless λ 1R cµ is large, the contributions proportional to the squared top Yukawa are the largest ones.
Assuming real couplings we get the following numerical expression in terms of LQ couplings:

A.11 LFV decays → γ, magnetic and electric dipole moments
We analyze in this section another set of LFV interactions, namely τ → µγ, τ → eγ and µ → eγ. Following [145], we can parametrize the α (p) → β (p − q)γ(q) vertex in terms of Lorentz invariant form factors: (A.108) which involves functions of p 2 , q 2 and pq. On-shell one has p 2 = m 2 α , q 2 = 0, pq = (m 2 α − m 2 β )/2, so that the form factors are just constants. In terms of the form factors the observables we are interested in can be written as where d and a stand for electric and (anomalous) magnetic dipole moment, respectively. In terms of LEFT coefficients we have In terms of parameters of the UV model, the leading terms are given by

JHEP01(2021)138
Particularly relevant for our model is the τ → µγ constraint. Assuming TeV-scale masses (to fix the logarithmic dependence on m 1 ) we get this approximate numerical expression: The experimental 95% CL bounds are [74] Br(µ → eγ) < 5.00 × 10 −13 , (A.115) The imaginary and real part of the diagonal terms T R αα define respectively the Electric Dipole Moment (EDM) d α and the anomalous magnetic moment a α for the three charged leptons, as seen in eq. (A.109). For the particularly relevant cases of the muon and electron magnetic moments, good numerical expressions, for TeV-scale masses and real couplings, are corresponding to a 3.7σ discrepancy: ∆a µ = (2.79 ± 0.76) × 10 −9 . Updated measurements will be available from the E989 experiment at Fermilab and E34 at J-PARC, which aim to reduce the experimental uncertainty by a factor of four.

JHEP01(2021)138
Other potentially interesting observables are quark EDM, which receive a strong bound via the constraint on the neutron EDM. Since these are not directly relevant for the parameter space relevant for our fits, we do not consider them in this work. A comprehensive analysis of EDMs in scalar leptoquark models is given in ref. [146]; for a recent analysis in the context of a (vector) U 1 -leptoquark model see ref. [147].

A.12 Z boson couplings
The couplings of the Z boson have been measured very precisely at LEP 1. At one-loop, the LQ model generates contributions to these very well measured quantities, which pose strong constraints on the model. The RGE-induced contributions in models aimed at addressing the B anomalies have first been studied in refs. [24][25][26]. Here we include the effect of finite corrections from the one-loop matching.
The pseudo-observables corresponding to the pole couplings of the Z boson to fermions correspond to these effective Lagrangian couplings: . The measurements of these pseudo-observables and the predictions for the SM contributions can be found in [99]. Also often used is the effective number of neutrino species at the Z peak, which depends on the couplings to neutrinos as (A.120) where g Z,SM ν ≈ 0.502. The latest update on the extraction of N ν from LEP data is given in [100]. We collect in table 3 the limits used in our fit.
Working at one-loop accuracy, there are two possible contributions to these pseudoobservables in our setup: tree-level contributions from the SMEFT operators, 13 one-loop matrix elements from the tree-level generated semileptonic operators. The result is: ij (µ) , Hq ] (1) ij (µ) , ij (µ) , (A.121) 13 There might be some indirect contributions via modifications to GF , but are negligible in our model.

JHEP01(2021)138
where µ is the scale at which the operators are evaluated, C X , and we included only one-loop matrix elements from up-type quark loops, since the top quark gives the dominant effect. The expressions for the one-loop matching of the SMEFT operators to the LQ model can be found in [58]. The loop functions are given by      Table 10. LEFT dipole (dimension five) operators.  -53 -