Palatable Leptoquark Scenarios for Lepton Flavor Violation in Exclusive $b\to s\ell_1\ell_2$ modes

We examine various scenarios that involve a light ${\cal O}(1 TeV)$ leptoquark state and select those which are compatible with the current experimental values for $\mathcal{B}(B_s\to \mu\mu)$, $\mathcal{B}(B\to K\mu\mu)_{\mathrm{large}-q^2}$, $R_K=\mathcal{B^\prime}(B\to K\mu\mu)/\mathcal{B^\prime}(B\to Kee)$, and which lead to predictions consistent with other experimental data. We show that two such scenarios are phenomenologically plausible, namely the one with a doublet of scalar leptoquarks of hypercharge $1/6$, and the one with a triplet of vector leptoquarks of hypercharge $2/3$. We also argue that a model with a singlet scalar leptoquark of hypercharge $1/3$ is not viable. Using the present experimental data as constraints, it is shown that the exclusive lepton flavor violating decays, $\mathcal{B}(B_s\to \mu\tau)$, $\mathcal{B}(B\to K \mu\tau)$ and $\mathcal{B}(B\to K^\ast \mu\tau)$, can be as large as $\mathcal{O}(10^{-5})$.


Introduction
The experimental searches of New Physics (NP) at the LHC are of great importance in our quest for solving the hierarchy and flavor problems. Rare decays and low-energy precision measurements play a crucial role in searching for the effects of physics beyond the Standard Model (BSM) because they are complementary to direct searches and even allow us to probe energy scales well above those that can be probed through them. Among the various types of rare decays, the lepton flavor violating (LFV) ones are particularly appealing because they are absent in the Standard Model (SM), and their detection would represent a clean signal of NP.
Theoretically, the LFV decays are predicted in various NP scenarios, such as in models with heavy sterile fermions [1], in models with a peculiar breaking of supersymmetry [2], in models involving an extra Z -boson [3], and in multiple scenarios founded on the paradigm of existence of one or more leptoquark (LQ) states [4]. Most of the current experimental efforts in searching for the effects of LFV are focused onto lepton decays, 1 → 2 γ, 1 → 3 2 (with 1,2 ∈ {e, µ, τ }), and on the µ−e conversion in nuclei [5]. The experimental sensitivity to probe these decays is expected to be improved in the years to come by several orders of magnitude. The LFV decays of hadrons are complementary to the purely leptonic decays because they represent a convenient ground to test LQ models, and they also represent a different environment to test the result of the above-mentioned leptonic modes. For example, in the LQ models, which will be discussed in this paper, the rates of the LFV decays of hadrons are much larger than in radiative and three-body lepton decays.
In this paper we will focus on the exclusive LFV modes, B s → 1 2 and B → K ( * ) 1 2 . At present, the most constraining experimental bound is B(B s → e ± µ ∓ ) exp < 1.1 × 10 −8 , which has been recently improved by an order of magnitude [6]. The only dedicated experimental searches with a τ -lepton in the final state have been made by the BaBar Collaboration which reported the upper limits B(B + → K + e ± τ ∓ ) exp < 3.0 × 10 −5 and B(B + → K + µ ± τ ∓ ) exp < 4.8 × 10 −5 [7]. The modes with one τ -lepton in the final state are phenomenologically more appealing because the relevant couplings are less constrained by direct experimental limits and for that reason, in the following, we will focus on B s → µτ and B → K ( * ) µτ .
One of the main motivations to study transitions based on b → s 1 2 comes from the recent observation of lepton flavor universality (LFU) violation in 1 where the squared dilepton mass is integrated in the bin q 2 ∈ [1, 6] GeV 2 . LHCb found [8] R exp K = 0.745 +0.090 −0.074 (stat) ± 0.036(syst), which is 2.4σ lower than the SM prediction R SM K = 1.00(1) [9,10]. Since the hadronic uncertainties largely cancel in the ratio, if confirmed, this result would be an unambiguous manifestation of NP. Another experimental observation of LFU violation has been made in the processes mediated by the charged-current interaction, namely, which turned out to be 2σ (for R D ) and 3.4σ (for R D * ) larger than predicted by the SM [11][12][13][14][15]. The possibility of connecting the observed LFU violation in neutral-and charged-current processes triggered many interesting theoretical works, mostly based on introducing various LQ states or new gauge bosons (Z and W ) [16,[19][20][21][22]. Remarkably, in most of the proposed scenarios, LFV in meson decays is induced at observable rates [17,18]. In a previous work, we computed the general expressions for the full angular distribution of B s → 1 2 and B → K ( * ) 1 2 decay modes, and explored a generic Z model, as well as the possibility of a Higgs induced LFV [24]. In this paper we make a phenomenological analysis of the exclusive b → s 1 2 decays in the LQ models which are compatible with the observed R K .
The remainder of this paper is organized as follows. In Sec. 2 we remind the reader of the effective approach to b → s 1 2 decays and briefly summarize the findings of Ref. [24]. In Sec. 3 we catalogue the plausible LQ models which are then scrutinized in Sec. 4 where we select the models compatible with the available b → sµµ data. The models selected in that way, are then subjected to a phenomenological analysis in order to obtain upper bounds on the LFV b → s 1 2 exclusive decay modes presented in Sec. 5. We summarize and conclude in Sec. 6. Most of the auxiliary formulas used in this paper are collected in the Appendices.

Effective Approach
The most general dimension-six effective Hamiltonian describing the LFV transitions b → s − 1 + 2 , with 1,2 ∈ {e, µ, τ }) is defined by [25] where C i (µ) and C 1 2 i (µ) are the Wilson coefficients, while the effective operators relevant to our study are defined by in addition to the electromagnetic penguin operator, O 7 = e/(4π) 2 m b (sσ µν P R b)F µν . The chirality flipped operators O i are obtained from O i by the replacement P L ↔ P R . Using the above Hamiltonian, one can easily compute the amplitudes for B s → − 1 + 2 and B → K ( * ) − 1 + 2 decays, and derive the expressions for the corresponding decay rates [24,26], also summarized in Appendix A.3 of the present paper. In the SM, the LFV Wilson coefficients, C 1 2 i , are zero due to lepton flavor conservation. The inclusion of neutrino masses can induce LFV at the loop-level, but with negligible rates due to the smallness of the neutrino masses. In the following we shall assume that the only source of LFV is NP and hence in order to compute the Wilson coefficients for b → s 1 2 ( 1 = 2 ) we need to specify the NP model.
To simplify our notation, in what follows, we will denote C i ≡ C 1 2 i ( 1 = 2 ) when confusion can be avoided. Notice, however, that: , which is particularly the case in the LQ models in which LFV occurs through tree-level diagrams; (ii) in some situations, even if C 1 2 i = C 2 1 i , ∀i, one can still generate an asymmetry between LFV modes with different lepton charges, e.g.
Before discussing the issue of lepton charge asymmetry one must first observe LFV, that is why we will here combine the two charged modes, namely, . We should also emphasize that the LFV channels we consider respect a hierarchy which depends on the adopted NP scenario. If the LFV is generated only by the (pseudo-)scalar operators, the lifted helicity suppression implies

(µ) in Leptoquark Models
Leptoquarks are colored states that can mediate interactions between quarks and leptons. Such particles can appear in Grand Unified Theories [27] and in models with composite Higgs states [28], as recently reviewed in [4]. In general, a LQ can be a scalar or a vector field which in turn can come as an SU (2) L -singlet, -doublet or -triplet [29]. In the following we assume that the SM is extended by only one of the LQ states and list the models which can be used to describe the b → s 1 2 processes. We adopt the notation of Ref. [30] and specify the LQ states by their quantum numbers with respect to the SM gauge group, (SU (3) c , SU (2) L ) Y , where the electric charge, Q = Y + T 3 , is the sum of hypercharge (Y ) and the third weak isospin component (T 3 ). Throughout this paper the flavor indices of LQ couplings will refer to the mass eigenstates of down-type quarks and charged leptons, unless stated otherwise. In other words, the left-handed doublets are defined as where V and U are the Cabibbo-Kobayashi-Maskawa (CKM) and the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrices, respectively, u L , d L , L are the fermion mass eigenstates, whereas ν L stand for the massless neutrino flavor eigenstates. Since the tiny neutrino masses play no role for the purpose of this paper, we may choose the PMNS matrix to be the unity matrix, U = 1. Right-handed field indices will always refer to the mass eigenbasis.

Scalar Leptoquarks
We first consider three scalar LQ scenarios. Two of them can modify the b → sµτ decay at tree-level without spoiling the proton stability. A third scenario, which might destabilise the proton and will also be discussed below, generates a contribution to b → sµτ only at the loop-level.
This model has already been used to provide a viable explanation of R exp K in Ref. [19]. The Yukawa Lagrangian reads, where ∆ ≡ iτ 2 ∆ * is the conjugated SU (2) L doublet, g L is a generic matrix of couplings.
In the second line we explicitly write the terms with ∆ (−1/3) and ∆ (2/3) where the superscripts refer to the electric charge of the Y = 1/6 LQ states. The masses of these two states can in principle be different, but for simplicity we will assume them to be equal, After integrating out the heavy fields, this model gives rise to the chirality flipped operators in the b → s − 1 + 2 effective Hamiltonian (4), and the corresponding Wilson coefficients are given by where v = 246 GeV is the well known electroweak vacuum expectation value The relevant Lagrangian for this model reads, where g R denotes the Yukawa coupling matrix, while again the superscript in ∆ (5/3) and ∆ (2/3) refers to the electric charge of the two mass degenerate LQ states. After integrating out the heavy fields, the only non-zero Wilson coefficients relevant to b → s − 1 + 2 are, corresponding to the chirality non-flipped operators in Eq. (4).
Being an electroweak singlet, this scalar LQ model is the simplest one. Its Lagrangian is given by, where the superscript C stands for the charge conjugation. 3 In addition to terms shown in (11) one could also write terms involving diquarks. To avoid conflict with the proton decay bounds those couplings must be negligible, so we set them here to zero. The Wilson coefficients contributing to b → s − 1 + 2 are generated at loop level and, to one-loop order, they are given by [21,33]: where f (x t ) = 1 + 3

Vector Leptoquarks
Vector LQs naturally appear in Grand Unified Theories as gauge bosons of the unified gauge group or in composite Higgs scenarios. Here we shall focus on relatively light vector LQs that can contribute to LFV b → s 1 2 transitions. There are three vector LQ states that could participate in d i → d j 1 2 processes via dimension-4 tree-level couplings to the SM fermions. These states, U 3 , V 2 , and U 1 , are respectively a weak triplet, doublet, and singlet under the SU (2) L . Here we do not consider the weak doublet vector V 2 , because its nonzero couplings might jeopardize the proton stability [30].
Being a SU (2) L singlet U 1 can couple to both the left-handed and right-handed fermions, namely, where x LL ij and x RR ij are the matrices of Yukawa couplings. The non-chiral nature of U 1 couplings to fermions induces both vector and scalar Wilson coefficients: where m U is the LQ mass. Obviously this model provides a rich ground for phenomenology because it gives rise to both the chirality flipped and non-flipped operators.

U 3 ≡ (3, 3) 2/3
The only possible dimension-4 gauge invariant interaction between U 3 , a weak triplet of LQ states, and fermions reads where by x LL ij we denote the couplings to the SM quarks and lepton doublets. The absence of diquark couplings guarantees that this LQ state is both baryon and lepton number conserving. The Wilson coefficients for the b → s 1 2 Hamiltonian are in the case of U 3 obtained by keeping nonzero only the x LL couplings in Eq. (15). For completeness we write them explicitly, where we assume the degeneracy among the charge eigenstates of the triplet, m U ≡ 4 Confronting LQ models with data 4.1 Which leptoquark model?
In this Section, from the list of models spelled out in Sec. 3, we select those that are consistent with the experimental data for the exclusive b → sµµ processes. To that end we require the consistency with the measured [36], since in that region of q 2 ∈ [15, 22] GeV 2 the hadronic uncertainties are well controlled by means of numerical simulations of QCD on the lattice [37]. Furthermore, by choosing q 2 ∈ [15, 22] GeV 2 we also avoid the narrow charmonium resonances so that we can rely on the quark-hadron duality [38]. This is also the reason why we use only B(B s → µ + µ − ) and B(B → Kµ + µ − ) high−q 2 to derive the constraints on possible NP contributions: (a) They involve the decay constant f Bs and a few form factors at large q 2 's, the quantities which have been accurately computed by means of lattice QCD; (b) These quantities are not plagued by uncertainties related to the cc-resonances. Theoretical description of B → K * µ + µ − decay entails more hadronic form factors and more assumptions are needed when confronting theory with experiment. Similarly, a recently observed discrepancy between theory and experiment regarding B(B s → φµ + µ − ) [39] relies on the assumption of validity of the light cone QCD sum rule estimates of the hadronic form factors at low and intermediate q 2 's for which the systematic uncertainties are hard to assess. Furthermore, for the purpose of our study, two measured quantities are sufficient to bound the NP couplings and we choose those which require the least number of assumptions, i.e. B(B s → µ + µ − ) and B(B → Kµ + µ − ) high−q 2 . Notice, however, that our results are compatible, at a (2 ÷ 3)σ level, with the global fit analyses made in Ref. [40].
• The ∆ (1/6) -model was studied in Ref. [19] where it has been shown that the required consistency with the measured B(B s → µ + µ − ) exp and B(B → Kµ + µ − ) exp high−q 2 results in the R K values fully compatible with experiment. 5 In this scenario (C 9 ) = −(C 10 ) , and from the combined fit of B(B s → µ + µ − ) and B(B → Kµ + µ − ) high−q 2 we obtain to 2σ accuracy, where the results in the first and second interval are obtained by using the hadronic form factors computed in lattice QCD in Ref. [42] and Ref. [43], respectively. Notice that the first interval coincides with what is obtained in Ref. [19]. In what follows we will choose 4 The quoted result was obtained after combining the CMS and LHCb results in Ref. [34]. Since then the ATLAS Collaboration also measured this same decay mode and reported B(B s → µ + µ − ) exp = (0.9 +1.1 −0.8 ) × 10 −9 [35]. Combining this result with the previous two would probably exacerbate the discrepancy with the SM but to do it properly one should combine the statistical samples of all three mentioned collaborations. 5 This model has also been considered in Ref. [41] to explain R exp K but by enhancing B(B → Kee), instead of decreasing B(B → Kµµ). which covers both of the above intervals.
• For the ∆ (7/6) -model, the NP Wilson coefficients satisfy the relation C 9 = C 10 . To scrutinize this scenario, we perform a fit of C 9,10 assuming them to be real and by applying the strategy described above. The allowed regions at 1σ and 2σ are shown in Fig. 1 where we also draw the lines corresponding to C 9 = ±C 10 . Clearly, the line C 9 = C 10 does not touch the region allowed by B(B s → µ + µ − ) exp and B(B → Kµ + µ − ) exp high−q 2 to 2σ. For that reason we discard the scenario ∆ (7/6) from further discussion. 6 Figure 1: Regions in the plane (C µµ 9 , C µµ 10 ) that are in agreement with the experimental values of B(B s → µµ) and B(B + → K + µµ) large q 2 at 1σ (dark green) and 2σ (light green) accuracy. The black point represents the SM prediction. The dashed lines correspond to the scenarios C µµ 9 = −C µµ 10 and C µµ 9 = C µµ 10 .
• The ∆ (1/3) -model generates both C 9 −C 10 and C 9 +C 10 . In other words, this model is compatible with exclusive b → sµµ data only if the combination with the plus sign is suppressed. One can achieve that by taking (g R ) ij ≈ 0 in Eq. (13), or by fine-tuning the terms in the same equation to get C 9 + C 10 ≈ 0. We will discuss more extensively this model in the next subsection.
• The U 3 -model was already studied in Ref. [20] where it has been shown that it can provide a simultaneous explanation of R exp K and R exp D ( * ) . In this scenario C 9 = −C 10 , which is in agreement with the measured B(B s → µ + µ − ) exp and B(B → Kµ + µ − ) exp high−q 2 , as shown in Fig. 1. To 2σ accuracy we find As before, we will use to cover the results obtained by using both sets of lattice QCD form factors.
• A similar discussion applies to the U 1 LQ in the limit of |x RR ij | |x LL ij | in Eq. (15). We should emphasize, however, that this model is not phenomenologically sound because the couplings to τ 's are very poorly constrained by data. Since one cannot compute the loop corrections with vector LQs without specifying the ultraviolet (UV) completion or providing an explicit UV cut-off, no constraints from B s −B s and τ → µγ can be used. While in the U 3 -model one can obtain a constraint by requiring the consistency between the tree-level process B(B → Kνν) and the upper experimental bound, in the U 1 -model the latter process is induced only at the loop-level and the experimental result is not constraining anymore. For this reason we prefer to discard the U 1 -model and focus only on the scenario with U 3 , as far as vector LQs are concerned.
In summary, of all the LQ models discussed in Sec. 3, we were able to discard in this section the scalar ∆ (7/6) and the vector U 1 . The remaining models can now be used to explore what are the allowed b → sµτ LFV decay rates. Before passing into that part, we revisit the ∆ (1/3) -model.

On (in)viability of
The interest in the ∆ (1/3) -model has been recently revived after the authors of Ref. [21] claimed that it can simultaneously accommodate the observed LFU violation in R K and R D ( * ) . In this section we revisit that claim and find that R K actually cannot be made significantly smaller than one, without running into serious difficulties with other measured quantities.
As we already mentioned above, the effective coefficients C 1 2 9 ± C 1 2 10 are loop induced and the results are given in Eqs. (12)- (13). To suppress the combination with the plus sign, which is disfavored by the current b → sµµ exclusive data, we set g R = 0 in Eq. (11). 7 The left-handed couplings are assumed to have the following flavor structure 7 Another possibility is to impose the cancellation between the terms in Eq. (13) but since the righthanded couplings lift the helicity suppression in leptonic decay rates one runs into problems with Γ(K → µν) and Γ(D (s) → µν), and can also significantly enhance Γ(τ → µγ), cf. Eq. (44) in Appendix A.1 of the present paper. In other words, the g R couplings are tightly constrained by data to be small, which justifies the approximation g R = 0.
where the first matrix connects down-type quarks to neutrinos, and the second up-type quarks to charged leptons. 8 Notice that the couplings to the first generation of leptons are assumed to be zero since they are tightly constrained by data. In our approach, the couplings to the d quark are also set to zero, since their nonzero values can induce significant contributions to processes such as K → πνν and the K −K mixing. Even though we fix (g L ) dµ = (g L ) dτ = 0, from Eq. (22), we see that the couplings to the u quark are generated via the CKM matrix. Therefore, constraints from the kaon and the D-meson sectors should be taken into account too. Besides those, one of the most important constraints comes from the B s −B s mixing amplitude. We obtain ∆m th where  [45,46]. A pecularity of this LQ model is the modification of the (semi-)leptonic decays of pseudoscalar mesons. We define the effective Lagrangian of the transitions d → u ν as where u and d stand for a generic up-and down-type quark flavor. Using this Lagrangian one can compute the (semi-)leptonic decay rates for the specific channels, e.g. B → D ν and B → ν , as described in Appendix A.1. After integrating out the LQ state, we obtain the effective coefficients which can be inserted into the expressions for decay rates, explicitly given in Appendix A.1, cf. Eqs. (46)- (47). We should stress that LQs can induce new contributions in which the neutrino has a different flavor from the charged lepton. One should therefore sum over the unobserved neutrino flavors in order to compare with the experimentally measured rates, e.g. B(B → D ν) = B(B → D ν ), with ∈ {µ, τ }.
Considering the ansatz given in Eq. (22) for the Yukawa matrix, the relevant leptonic modes for our study are K → µν, D s → (µ, τ )ν, and B → τ ν. We consider the experimental values given in Ref. [47] and we use the values for the decay constants computed in lattice QCD, summarized in Appendix A.4. Other observables, such as B(τ → µγ) and B(D 0 → µµ), for which the formulas are given in Appendix A.1, give important constraints as well.
In addition to the above-mentioned constraints, compactly collected in Tab. 1, we impose the conditions of perturbativity, |(g L ) ij | < 4π, and look for the points which would simultaneously satisfy the observed R exp K and R exp D . As it can be seen in Fig. 2 (left panel), we were not able to find couplings (g L ) ij that would result in values for R D consistent with experiment, R exp D = 0.41(5) [11,12]. In fact, the couplings of ∆ (1/3) to the muon, which are necessary to get R K < 1, are large enough and push R D to values smaller than the SM one. 9 Furthermore, even though we were able to find points that give acceptable values for R K , we find that the selected points are in conflict with as depicted in the right panel of Fig. 2. Although R µ/e D has not been experimentally established, values of R µ/e D ≈ 1.05 seem already implausible, since the B(B → Deν) and B(B → Dµν) data have been successfully combined in B-factory experiments to extract G(1)|V cb |. In Ref. [22] it was even argued that such a deviation from lepton flavor universality cannot be larger than 2%. In Fig. 2, however, we see that the points selected to satisfy the observed R exp K result in R µ/e D > 1.8, a large departure from one. We, therefore, conclude that one cannot accommodate the experimental value R exp K without producing a huge enhancement of B(B → Dµν), which implies R D R SM D and unacceptably large R µ/e D . 10 We were insisting on the agreement with the experimental value R exp K . If, instead, one wants to accommodate the experimental value R exp D with this model [44], the resulting values of the couplings to muon are either zero or far too small to explain R exp K . There are several sources of disagreement between the present paper and Ref. [21], which we comment on in detail. First of all, their best fit values for R D ( * ) were obtained in a setup where NP only couples to τ 's. This assumption is not justified, since a putative explanation of R K through loops necessarily requires large couplings to muons, which then modifies the denominator in R D , and produces excessively large R µ/e D , as discussed above. Secondly, the conditions for the LQ parameters given by Eqs. (12) and (17) of Ref. [21] are necessary but not sufficient to reproduce results consistent with B(B → Kνν) exp and R exp K , respectively. Finally, the bounds from processes involving the up-quark, such as K → µν, D s → (µ, τ )ν and B → τ ν, have been tacitly neglected in their work. As mentioned before, one cannot completely avoid the first quark generation, see Eq. (22). 9 In obtaining R D we used the B → D form factors recently computed in lattice QCD [48], which also give R SM D = 0.286 (12). 10 In obtaining R

Results and discussion
In this section we present our predictions for the remaining viable scenarios, namely the LQ models ∆ (1/6) and U 3 .

∆ (1/6) ≡ (3, 2) 1/6
In our analysis of the ∆ (1/6) -model we assume the LQ couplings to the first generation of fermions to be negligible because they are most tightly constrained by the experimental limits on µ − e conversion on nuclei, on atomic parity violation, on B(K L → µe), and on B(B s → µe). Therefore, the matrix of couplings (g L ) ij in Eq. (7) is assumed to have the form where the entries are considered to be real. The product (g L ) sµ (g L ) * bµ = 0 is then fixed by Eqs. (8) and (19). To constrain the couplings to τ 's we use the information from B s −B s mixing amplitude, and the experimental limits on B(B → Kνν) and on B(τ → µφ), see also Tab. 1. To that end we recall the expression for the mass difference of the B s −B s system with respect to the SM one, ∆m th Bs ∆m SM
Other observables, such as leptonic decays of φ(nS) and Υ(nS) mesons, could in principle provide useful bonds on the LQ couplings but the derived limits turn out to be much less significant at this point. Flavor conserving leptonic decays, such as Υ(nS) → τ τ , are dominated by the very large tree-level electromagnetic contribution which undermines their sensitivity to NP. LFV modes, such as Υ → µτ , are in principle sensitive to the LQ contribution but the current experimental limit B(Υ(1S) → µτ ) exp < 6 × 10 −6 [54] is still too weak to be useful. We have also checked that in this model the upper experimental limit on B(τ → µγ) does not provide any additional constraint because of the cancellation of terms ∝ 1/m 2 ∆ in the loop function, see Eq. (42) in Appendix A.1. Using the constraints discussed above, and summarized in Tab. 1, in addition to the perturbativity prerequisite, which we here decide to be |(g L ) ij | ≤ 1, we were able to scan the parameter space and find points which are consistent with all the requirements. With the points selected in that way, we then compute the Wilson coefficients (C µτ 9 ) = −(C µτ 10 ) by Quantity m ∆ = 1 TeV m ∆ = 5 TeV m ∆ = 10 TeV B(B s → µτ ) < 1.0 × 10 −5 < 3.0 × 10 −6 < 1.8 × 10 −7 B(B → Kµτ ) < 1.1 × 10 −5 < 3.4 × 10 −6 < 2.0 × 10 −7 B(B → K * µτ ) < 2.0 × 10 −5 < 6.1 × 10 −6 < 3.7 × 10 −7 where in evaluating the above ratios we used the central values of the hadronic parameters.
Since most flavor observables are studied at tree-level, our predictions depend only on the ratios (g L ) d /m ∆ with d = d, s, b and = µ, τ . The only exception is ∆m Bs , which exhibits a mild dependence on m ∆ . We see from the results presented in Tab. 2 that the bounds become tighter for larger m ∆ which is a consequence of the condition |(g L ) ij | ≤ 1.  [46] and B(B + → K + µ ± τ ∓ ) exp < 4.8 × 10 −5 [7].
We next check on the correlation between B(B → Kµτ ) and R νν for a benchmark m ∆ = 1 TeV. We see from Fig. 3 that the compatibility of R νν with the SM value, R SM νν = 1, would not exclude the possibility of rather large branching fractions for the LFV decay modes. However, if the value of R νν is found to be significantly below or above 1 then this analysis provides both an upper and a lower bound for B(B → Kτ µ). For example, if Furthermore, we see that the current experimental bound, R exp νν ≤ 4.3, is an efficient constraint but it ceases to be so if we consider m ∆ 3 TeV. Finally, it is interesting to note that the BaBar bound, B(B → Kµτ ) exp < 4.8 × 10 −5 , is only an order of magnitude larger than the upper bound we obtained in our analysis. We hope that the high statistics LHCb data, and especially those of the future Belle II, will be used to improve this as well as other bounds on similar decay modes. We reiterate that a bound on any of the modes considered above would be equally useful, since they are all related to each other through the ratios given in Eq. (33). Before closing this part of our analysis, we should mention that this LQ model has been extended to accommodate R exp D in Ref. [55], without changing the upper bounds on the LFV rates discussed in this paper but providing the lower ones which are O(10 −10 ).

Comment on B s → τ τ
Since the LHCb Collaboration is actively searching for the B s → τ τ events, we also show the correlation between B(B → Kµτ ) and B(B s → τ τ ) in Fig. 4. We see again that even if B(B s → τ τ ) is found to be consistent with the SM value, B(B s → τ τ ) SM = 7.7(5) × 10 −7 , one can still have a significant B(B → Kµτ ). This can be understood from Eq. (8), since B s → τ τ will be modified by NP only if the product (g L ) sτ (g L ) * bτ is nonzero, and therefore one should have both (g L ) sτ = 0 and (g L ) bτ = 0 to produce a deviation from the SM. On the other hand, the product (g L ) sµ (g L ) * bµ is already fixed to a nonzero value by b → sµµ exclusive data. Therefore, it is enough to have one of the couplings (g L ) sτ or (g L ) bτ different from zero in order to induce LFV in B (s) decays through the effective coefficients (C τ µ 10 ) or (C µτ 10 ) , respectively. For that reason, the LFV modes are more robust probes of the LQ couplings. Notice, however, that if LHCb can establish an upper bound on B(B s → τ τ ), this would certainly be more important than the other limits that we have used here to constraint the third generation couplings.
If, instead, the measured B(B s → τ τ ) turns out to be larger or smaller than B(B s → τ τ ) SM this model offers also a lower bound on B(B → Kµτ ), which is quite remarkable because in this way one can check the validity of this scenario. Moreover, we see also that if B(B s → τ τ ) 2 × 10 −5 is measured then this model not only cannot generate LFV, but it will be ruled out altogether.

U 3 ≡ (3, 3) 2/3
In a UV-complete renormalizable model containing U 3 at the lower end of the spectrum (and possibly accompanying light degrees of freedom), we expect the 3 × 3 matrix of LQ couplings x LL in Eq. (16) to be unitary. Indeed, if U 3 was a remnant of a larger multiplet of gauge bosons, one could consider a theory in the unitary gauge in which the UV divergences appearing in the flavour changing neutral processes are removed by the sums over the couplings, i.e. by the GIM mechanism, and for this technique to work we would require x LL to be unitary. 11 However, in the unitary case the couplings to e and/or first generation of quarks cannot be avoided and then the very strong null constraints from LFV searches apply. In this particular case, we have found that the upper bounds on µ − e conversion on Au nucleus and the one on B(K L → µe) already discard U 3 with unitary x LL as a viable scenario for explaining R exp K . Thus we continue along with our analysis using the non-unitary ansatz, and we completely avoid couplings to electrons by paying the price of being unable to unambiguously predict any loop-induced process. The ansatz below couples only the second and third generation charged leptons with the down-type quarks. However, as before, the CKM rotation will induce the couplings of u, c, t to the charged leptons and neutrinos, viz.
Observables that constrain this LQ have already been discussed in Ref. [20]. Lefthanded couplings of U alone allow for a single combination of the Wilson coefficients: As discussed in Sec. 4.1, by using the measured B(B s → µ + µ − ) exp and B(B → Kµ + µ − ) exp high−q 2 , we select the values of C µµ 9 = −C µµ 10 , to 2σ accuracy and obtain C µµ 9 ∈ (−0.76, −0.04), cf. Eq. (21). It is interesting to note that B(B s → µ + µ − ) exp and R exp K impose almost degenerate constraints on C µµ 9 . For example, the result for C µµ 9 = −C µµ 10 to 1σ, results in R K ∈ (0.74, 0.88), which fully agrees with R exp K . 11 See also the discussion in Ref. [56].
The three charge components of U 3 , with the same underlying couplings, allow for flavor effects both in the down/up-quarks and the charged lepton/neutrino sectors. Interestingly, the 2/3-charge component plays an important role in the charged current semileptonic processes and thus can be used to address the hint of LFU violation in R D ( * ) . Using the effective Lagrangian parametrization of the charged current (24) the left-handed current modification reads Through g V the U 3 vector LQ state modifies the overall normalization of the SM spectrum of B → D ν without modifying its shape, cf. Eq. (47) in Appendix A. Thus in the U 3 -model we vary the couplings within |x LL ij | < 4π and apply the abovementioned constraints, namely the one on C µµ 9 , and the one arising from R exp D , both at 2σ. 12 We fix the mass at m U = 1 TeV. Further constraints are coming from the purely thirdgeneration charged current and the up-type neutral current processes which are mediated by the U (5/3) 3µ eigencharge component. Among the charged current processes, the most efficient constraint comes from the measured B(t → bτ ν) exp = 0.096(28) [63]. The effect of U 3 here is an overall rescaling of t → bτ ν via coupling g V | b→tτ ν , cf. Eq. (37). We apply this constraint by eliminating any point in the parameter space that results in B(t → bτ ν) out of the 1σ region of B(t → bτ ν) exp . Once the above constraints are applied it appears that the constraints coming from the up-quark rare process, B(D 0 → µµ) exp < 6.2 × 10 −9 [47] (cf. Eq. (53)), and from B(τ → µφ) are redundant. Furthermore, we have checked that the LFU ratios, Γ(τ → Kν)/Γ(K → µν) and Γ(K → µν)/Γ(K → eν), are consistent with experiment.
Like in the previous subsection, with the couplings selected to verify the above constraints, we can compute various quantities. In particular, in Fig. 5 we show the correlation between R νν and B(B s → τ τ ), and again observe that the current experimental limit on R νν provides a stringent constraint on the parameter space. Note that B → Kνν is governed by the Q = −1/3 component which connects the down-type quarks to neutrinos.
Finally, on the right panel of Fig. 5 we show B(B s → µτ ) with respect to R νν , and we see that this model allows for large LFV decay branching fractions even if R νν and B(B s → τ τ ) are close to their SM values. We obtain for m U 1 TeV. The bounds on the other two LFV modes can be obtained by simply using Eq. (33), i.e.
The strong correlation between LFV and bi-neutrino modes is a signature of this model, should positive LFV signal or R νν > 1 emerge. On the other hand, for the B s → τ τ mode the left-hand panel in Fig. 5 demonstrates that it can be larger by about an order of 12 To explain R exp D in this model, |x LL bµ | should be 1, which is why we opted for |x LL ij | < 4π, as often used in the literature as a perturbativity requirement for the couplings. magnitude than its SM prediction. For R νν > 1 there is a possibility of strong suppression of B(B s → τ τ ). Finally, the observation of R νν significantly smaller than 1 is not compatible with the presence of U 3 as an explanation for R exp K and R exp D . We should emphasise that all the observables considered in this section have identical scaling with U 3 couplings and mass, thus the correlations presented in Fig. 5 are kept intact for different choices of m U , with the exception of perturbativity constraints which are becoming more stringent for increasing m U .

Summary
In this paper we discuss exclusive b → s 1 2 LFV decay modes in various models in which the existence of an O(1 TeV) leptoquark state is assumed. In doing so we use the effective field theory approach and derive the expressions for the Wilson coefficients in each of the considered scenarios, which then allows us to test them against the available experimental data for B(B s → µ + µ − ) exp and B(B → Kµ + µ − ) exp at large q 2 , and select those which lead to predictions of other processes that are consistent with current experimental data. These two quantities are used mostly because the theoretical uncertainties are under better control since the hadronic matrix elements have been estimated by means of numerical simulations of QCD on the lattice, and also because the high-q 2 region is above the very narrow ccresonances, so that the contribution of cc can be treated by using the quark-hadron duality.
We find that two models (scenarios), compatible with the above-mentioned experimental constraints, can accommodate the experimental hint on LFU violation, namely R exp K = B (B + → K + µµ)/B (B + → K + ee) < R SM K , and that they lead to the results in other processes that are compatible with experiment (and with the measured ∆m exp Bs , in particular). Moreover, we find that the current experimental bound R exp νν is particularly useful in cutting a significant fraction of the parameter space. The two scenarios that we select as plausible are: • Scenario with ∆ (1/6) , a doublet of (mass degenerate) scalar leptoquarks of hypercharge Y = 1/6, for which the NP contribution to B(B s → µ + µ − ) and B(B → Kµ + µ − ) comes from the left-handed couplings, and for which the bound B(τ → φµ) exp appears to be useful. Interestingly, we find that the LFV mode B(B → Kµτ ) can go up to O(10 −6 ) even if R νν = 1 and B(B s → τ τ ) exp = B(B s → τ τ ) SM , but that it is strictly different from zero for R νν = 1 and B(B s → τ τ ) exp = B(B s → τ τ ) SM . For example, for R νν ≈ 2, in this model we get B(B → Kµτ ) ∈ (0.2, 50) × 10 −7 . Notice that this model is experimentally verifiable: (i) it leads to R K * > 1 [19], and (ii) it does not allow R νν 0.6, nor B(B s → τ τ ) 2 × 10 −5 .
• Scenario with U 3 , a triplet of (mass degenerate) vector leptoquarks of hypercharge Y = 2/3, for which the NP contribution to B(B s → µ + µ − ) and B(B → Kµ + µ − ) comes from the left handed couplings, and for which the bound on R exp νν as well as B(t → bτ ν) exp provide the crucial constraints on its Yukawa couplings. Like in the previous case, the branching ratio B(B → Kµτ ) can be between zero and O(10 −6 ) even if R νν = 1 and B(B s → τ τ ) exp = B(B s → τ τ ) SM , but it is stricly non-zero for R νν > 1. This model too is experimentally verifiable: (i) it leads to R K * = 0.86 (12), thus R K * < 1 [41], and (ii) it can be ruled out if R νν < 1 and/or B(B s → τ τ ) 10 −5 .
As for the branching fractions of the similar LFV decay modes, they are related to B(B → Kµτ ) via Eq. (33), as discussed in Ref. [24]. Notice that we also examined the model with a singlet scalar leptoquark of hypercharge Y = −1/3 and found that it is not phenomenologically viable, i.e. that it cannot accommodate R K < 1 without running into serious phenomenological difficulties with other measured processes.

A Formulas and hadronic quantities
In this Appendix we collect the expressions and numerical values used in our analyses.
• B → τ ν and B → Dτ ν: Using the effective Lagrangian defined in Eq. (24) one can compute the leptonic and semileptonic decay rates for the specific channels, e.g. for B → D ν and B → ν . We obtain and dΓ where we emphasize that the NP contribution can allow for a different neutrino flavor = in the final state. The B → D form factors and the f B decay constant are defined in the Appendix A.3, and the phase-space functions c i (q 2 ) are given by [15] c + (q 2 ) = λ 3/2 ( q 2 , m B , m D ) 1 − 3 2 A.
• D 0 → µ + µ − : The branching ratio of the rare decay D 0 → µµ is sensitive to the CKM-rotated matrix x LL : • τ → µφ: The expression relevant to this LFV decay mode in the U 3 -model reads,

A.3 Exclusive b → s decay rates
We recall in this section the full decay rate expressions for the modes B s → 1 2 and B → K ( * ) 1 2 [24].
The branching fraction for the mode B s → 1 2 reads where the B s -meson decay constant f Bs is defined by 0 b γ µ γ 5 s B s (p) = ip µ f Bs .

A.4 Hadronic quantities
In this section we summarized the hadronic parameters considered in our analysis. The decay constants obtained by simulations of QCD on the Lattice are summarized in Tab. 3 [37]. In order to perform the fit of B(B → Kµ + µ − ) high−q 2 , we combine the B → K form factors which were precisely determined in the low-recoil region by Ref. [42], and more recently by Ref. [43]. Similarly, to compute R D we have used the B → D form factors recently computed in Ref. [48]. Finally, the LFV semileptonic decays B → K ( * ) 1 2 were analysed by employing the form factors of Ref. [61].