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

We perform a complete study of the low-energy phenomenology of $S_1$ and $S_3$ lepto-quarks, 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.

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 the leading radiative effect [24][25][26][27], however since the logarithm is often just of O(1), finite contributions can have a relevant impact.
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 a 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. Ref. [56] for a recent review of collider searches for S 1 and S 3 .
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 Sec. 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. The results for all scenarios considered are collected in Sec. 3 and a discussion on future prospects can be found in Sec. 4. We conclude in Sec. 5. In App. 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 2 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 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: [C (1) lq ] 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: 1. The one-loop matching for the S 1,3 model into the SMEFT, up to dimension-six operators, resulting by integrating out the two scalar leptoquarks at a scale of the order of their masses µ M ∼ M 1 , M 3 . The complete set of matching conditions, obtained with MS renormalization scheme, has been provided in [58]. 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 i , only the tree-level matching conditions from SMEFT to LEFT, and tree-level 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.  [77] (11.0 ± 4.0) × 10 −11 [78]

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 Table 1, we show the list of low-energy observables that we analyze, together with their SM predictions and experimental bounds. In App. 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

SM prediction
Experimental bounds ∆F = 2 processes App. A.7  [23] (279 ± 76) × 10 −11 [22,23] 2.9963 ± 0.0074 [95]  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 App. 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. For each scenario we get the CL regions in the plane of two real couplings, by profiling the likelihood over all the other couplings. We are often also interested in the values of some observables corresponding to these CL regions. To obtain this, we perform a numerical scan over all the parameter space 4 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.

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: (3.1) 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 described in the previous Section, with all the observables. For instance, if the couplings to muons are set to zero, neutral-current B-anomalies and the muon anomalous magnetic moment are automatically frozen to the corresponding Standard Model values and do not impact the final fit.
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.
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): Table 4: Summary of leptoquark models considered in this work. The third columns lists the couplings we allow to be different from zero in our global fit. The last three columns indicate whether the models provide a satisfying fit of each set of anomalies, respectively.
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.1 Single-leptoquark S 1

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: R(D), R(D * ), Br(B + c → τ + ν) (App. A.2), δg Z τ R (App. A.12), and |g τ /g µ | (App. A.8). The results from the fit, assuming real couplings, can be seen in Fig. 1. All observables in this case scale with λ/M , so we show the limits normalizing the leptoquark mass to 1 TeV. The left panel shows confidence regions for the two couplings. Notice that, even though the best-fit value for λ 1R cτ (indicated as a black dot) may shed doubts on perturbativity, a large region of smaller values is allowed by the fit. The dashed lines are 95% CL constraints from single observables, and help illustrate the role of each observable within the global fit.
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].
The model is successful in fitting the deviation in R(D ( * ) ) within the 1σ level, with smaller values of R(D * ) (or larger R(D)) preferred by the fit. Also the fit of ∆g Z τ R is somewhat improved with respect to the SM, which is in mild tension with the experimental value (see Table 3). Improved measurements of R(D ( * ) ) can test this setup due to the precise linear relationship among the two modes predicted by the model.
bµ . Right: 95% CL limits from individual observables in the plane of λ 1L bµ − λ 1L sµ , fixing M 1 = 6 TeV. In both plots the green region is the 1σ favourite one from ∆C µ 9 .
The best-fit point 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), c.f. Eqs. (A.4,A.7). B s -mixing, c.f. Eq. (A.72), and B → K ( * ) νν, c.f. Eq. (A.47), 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 Fig. 2 (left). For M 1 3 TeV it is possible to avoid these limits while having couplings still in the perturbative range (see also [38]). Nevertheless, even while marginally evading the B s -mixing constraint, the ∆C µ 9 deviation remains in ≈ 2σ tension with the bound on λ 1L bµ arising from the LFU limit from τ decays, |g µ /g e | (see App. A.8). The situation is illustrated in the right panel of Fig. 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: Thanks to the m t enhancement of the left-right contribution, c.f. Eq. (A.107), 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.112)): 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. [96]. More recently, in the updated versions of [97,98] 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 charm. 5 In our formalism this can be achieved aligning the S 1 couplings as 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 S 1 . The relevant couplings are 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 Banomaly), and λ 1R tµ , λ 1L bµ on the other hand (necessary to fit a µ ), sizeable contributions to τ → µγ are generated at one-loop, see Eq. (A.109). 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.01, corresponding to a tuning of approximately one part in 15. 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 Fig. 3. In the upper two panels we show the fit in the (λ 1L bµ , λ 1R tµ ) and (λ 1L bµ , λ 1R tτ ) planes, including We show the preferred regions in the planes of two couplings, where those not shown are profiled over, and the individual 2σ limits from the most relevant observables for illustrative purposes. The black dot corresponds to the best-fit point.
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 fully reproduced, and that τ → µγ requires a non-zero value of λ 1R tτ , as discussed above. In the lower panel we show the fit in the couplings contributing to R(D ( * ) ), (λ 1L bτ , λ 1R cτ ). Comparing to the allowed region in the same plane in the model studied in Fig. 1, we see that τ → µγ strongly reduced the allowed region, while still not preventing a very good fit of the charged-current B-anomalies.
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, with M 1 = 1 TeV, for Figure 4: Result from the fit in the S 3 model. In the upper panels we show the preferred regions in the planes of two couplings, where the two not shown are profiled over, and the individual 2σ limits from the most relevant observables for illustrative purposes. In the lower panels we show where the preferred region is mapped in the planes of the neutral and charged-current anomalies.

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 Fig. 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 chargedand 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.
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.27).
The most salient features of the fit are summarized in Fig. 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 Fig. 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.41). 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.4 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 Fig. 6.
In the first row of Fig. 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 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  Observable Measurement . 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 Fig. 6  (see Fig. 4 top-left). We conclude that all the anomalies in R(D ( * ) ), b → sµµ, and (g − 2) µ , can be completely addressed in this model, for perturbative couplings and TeV-scale leptoquark masses. The best-fit point, for M 1 = M 3 = 1 TeV, is found for

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 [103]: S,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 [104] and our one-loop matching of the SMEFT to the LQ model, we obtain (see App. 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 Figure 7: Limits on LQ potential couplings from oblique corrections and Higgs measurements.
In each panel, the other two couplings have been marginalised. The black point represents the best-fit point while the dashed blue contours are the prospects for 95%CL limits after HL-LHC.
negligible given the present experimental precision. The constraints on S and T from [99] are reported in Table 5. The contribution to the T parameter from the λ H13 coupling has 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 App. 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 Fig. 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 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 [105]. 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. 6 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 , c.f. 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 App. A.6); • B decays induced at parton level by b → sτ τ (see App. 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 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.

Observable Present limit
Belle II (5)50ab −1 LHCb Up.-II b → sτ µ observables Br(B + → K + τ ± µ ∓ ) < 3.3(5.4) × 10 −5 [80,81] 3.9 × 10 −6 O(10 −6 )  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 [106] and LHCb [107] are summarised in Table 6.
In Fig. 8 we show how the preferred parameter-space regions for the models S 1 +S (right) map in the plane of the branching fractions of the LFV decay B + → K + τ µ and the decay B + → K + τ τ (normalised to the SM value). 7 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 + S 3 (all) model. For the leptonic decay B s → τ + τ − at Belle II only the prospect at luminosity of 5ab −1 is available; it is not shown in the plots since it is weaker with respect to the semileptonic decays and correspond to a horizontal line at ∼ 1250. On the other hand, for the Upgrade II of LHCb, the prospected bound on Br(B s → τ + τ − ) (blue horizontal lines) is stronger and leads to constraints similar to the ones that we obtain from the B + → K + τ + τ − decay measured at 5ab −1 Belle II. In order to evaluate the constraining power of future Br(B s → τ µ) measurements, in Fig. 8, one could keep in mind that in our model we have Br(B s → τ µ)/Br(B + → K + τ µ) ≈ 0.89, at tree-level.

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 Bmeson 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 Fig. 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 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 Sec. 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 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 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.
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.

A Analysis of observables and pseudo-observables
We collect in this Appendix the details for all observables considered in this study. Whenever feasible and relevant, we emphasize the model independent steps of the calculations. In the derivations below, we employ several results available from the literature. Of particular relevance to us are the following references: • the complete one-loop matching equations of the S 1 + S 3 model onto the SMEFT [58] 8 , • the SMEFT one-loop renormalization group equations [61][62][63], • the one-loop matching equations of SMEFT onto LEFT [64], • the LEFT one-loop renormalization group equations [65].
We consider up to dimension-six operators in the effective theories. Operator bases are given in Refs. [59] and [110], for SMEFT and LEFT respectively. We give, for convenience, our numerical expressions in terms of m 1,3 ≡ M 1,3 /(1 TeV).

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 8 Notice that our previous work [58] employed the opposite sign convention for gauge couplings, with respect the present paper (cf. Eq. (2.2)). We switch our convention for consistency with Refs. [61][62][63], [64] and [65]. We also note that Refs. [61][62][63] employ a different convention for SM Yukawa couplings (i.e. y [61][62][63] = y † this work .) The expressions for the Wilson coefficients of the above operators in terms of the LEFT ones are Experimental measurements of the observables taken into account set constraints on the C sb 9(10) coefficients of O ( ) 9(10) . The only tree-level contributions come from [O V,LL ed ] sb . One-loop corrections might be important in a precision analysis, in particular in regions 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.3) are all vanishing at first approximation, apart from C sb 9(10) which can be obtained, using Eq. (A.4), as The relevant observables for these operators are the LFU ratios R K [12,14,111] and R K * [13,15], the branching ratio of B s → µ + µ − [19,20,112], the angular observable P 5 in B → K * µ + µ − [21], and branching ratios in other b → sµµ transitions [18,[113][114][115]. 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 lefthanded combination ∆C sbµµ with a correlation of ρ ≈ −0.5. 9 The numerical expressions we get in our setup are, assuming real LQ couplings: 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 [116] 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  [117] (while L V,LL νedu is not affected by this running), and by electroweak and Yukawa RGE from the electroweak to the leptoquark scale. In our 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:   [58]. Considering a NP scale around ∼TeV, an approximate expression of the Wilson coefficients at m b scale as a function of the LQ couplings, is In terms of low-energy coefficients at the m b scale, we obtain the relation C S L ≈ −7.  [118]: In the numerical analysis we use the complete expressions from [118], including quadratic terms. The global avarage of R(D) and R(D * ) [1][2][3][4][5][6][7][8][9][10][11] is taken from [67], the measurement of P τ (D * ) from [7], F D * L from [69], and Br(B + c → τ + ν) from [70]. They are summarised in Table 1 together with the SM predictions.
Approximate numerical expressions for R(D), R(D * ), and Br(B + c → τ + ν) assuming real LQ couplings: 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 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. [119]. 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.11) (but with = µ, e instead of the tau) is given by [120] 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.27) 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 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] R These ratios can be expressed in terms of 2 real parameters > 0 and η ∈ [−1/2, 1/2] [75]: (A.45) 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.43) 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: Let us now discuss Kaon decays. In the Standard Model, the coefficients of this four-fermion interaction are [77] [C SM,ds with X t = 1.48, X e c = X µ c = 1.053 × 10 −3 , and X τ c = 0.711 × 10 −3 . Since C R ≈ 0 in our model, the branching ratios for the K + → π + νν and K L → π 0 νν decays can be expressed as where Br(K + → π + ν e ν e ) SM = 3.06 × 10 −11 , Br(K + → π + ν τ ν τ ) SM = 2.52 × 10 −11 , Br(K L → π 0 νν) SM = 3.4 × 10 −11 . (A.51) The experimental bounds are [74,78] 10 Br(K + → π + νν) = (11.0 +4.0 −3.5 ) × 10 −11 , Br(K L → π 0 νν) < 3.57 × 10 −9 (95%CL) . (A.52) Assuming that the contribution in ν τ is the dominant one we get the following approximate numerical expression:  10 NA62 results on K + → π + νν have been updated in a presentation at the ICHEP2020 conference.
We consider the leptonic branching ratios: 54) and the semileptonic ones: At partonic level, these processes are induced by the same low-energy operators as b → s , see Eqs. (A.3,A.4), 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.3,A.4), C (9−10) ≡ C 9 − C 10 and we have Numerically, we get: (A.61) The present 95% CL limit on this branching ratio is 6.8 × 10 −3 [109], 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 [123] 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 [123]. The leading contribution is given by Present limits and future prospect for the observables discussed in this Section are reported in Table 6.
They receive NP corrections, at one-loop level in the S 1,3 model, from the four-quark SMEFT operators O In the S 1 + S 3 model, the non-vanishing coefficients are: 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 [82]. 11 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 combina-tions are We remark that the bounds in Ref. [83] 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 [84] 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.73) 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.9 LFV decay τ → µφ The relevant effective operators for the τ → µss transition in the LEFT are: [L S,RL ed ] τ µss (τ R µ L )(s R s L ) + h.c.

(A.81)
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 QED RG contribution) and keeping only couplings to third generation (as well as λ 1L sµ , which can be large in some scenarios). In this case we get: where terms with coefficients smaller than 10 −4 are not shown. Among the several observables testing the τ → µss contact interaction, the branching ratio of τ → µφ gives the most stringent constraints [124]. It is given by (see also [45,124]) [45,74], and, at tree-level, 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. 91) A.10 LFV decay τ → 3µ In the low-energy EFT the three body LFV decay τ → 3µ is induced by four-lepton and lepton dipole operators: where we explicitly indicate only the ones that are generated at dimension-six level in the SMEFT. In our LQ model, all these operators arise at loop level. 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.
The branching ratio for the three body LFV decay τ → 3µ can be expressed as [125] Br(τ → 3µ) = m 5 Above, X γ parametrizes the contribution from the dipole operator O eγ and is given by ] µτ µµ holds. Assuming real couplings we get the following numerical expression in terms of LQ couplings: If the S 1 couplings to RH fermions vanish, the contributions to the process are much smaller, the leading terms being: Presently, the experimental bound on this LFV decay, at 95% C.L. is [74] Br(τ → 3µ) < 2.5 × 10 −8 .
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. [126]; for a recent analysis in the context of a (vector) U 1 -leptoquark model see Ref. [127].

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 Ref. [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 [94]. Also often used is the effective number of neutrino species at the Z peak, which depends on the couplings to neutrinos as (A. 116) where g Z,SM ν ≈ 0.502. The latest update on the extraction of N ν from LEP data is given in [95]. 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 12 , one-loop matrix elements from the tree-level generated semileptonic operators. The result is: ij (µ) , Hq ] (1) ij (µ) , ij (µ) , 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 12 There might be some indirect contributions via modifications to G F , but are negligible in our model.

A.13 Oblique corrections and Higgs couplings
Let us now consider a set of constraints for the potential couplings between leptoquark and the SM Higgs. Since these are universal corrections we focus on constraints independent on the LQ couplings to fermions. Considering that all effects related to the couplings in the potential arise at one-loop, and that we also have a scale suppression from the TeV-scale LQ masses, it is clear we need to consider only high precision constraints for universal theories: EW precision observables, and Higgs processes which arise at one-loop in the SM.