The fate of $V_1$ vector leptoquarks: the impact of future flavour data

Motivated by the recent experimental progress on the $B$-meson decay anomalies (in particular the angular observables in $B\to K^\ast\mu\mu$), we rely on a simplified-model approach to study the prospects of vector leptoquarks in what concerns numerous flavour observables, identifying several promising decay modes which would allow to (indirectly) probe such an extension. Our findings suggest that the confirmation of the $B$-meson decay anomalies, in parallel with positive signals (at Belle II or LHCb) for $\tau\to \phi \mu$, $B_{(s)}$-meson decays to $\tau^+ \tau^-$ and $\tau^+ \mu^-$ ($\tau^+ e^-$) final states, as well as an observation of certain charged lepton flavour violation decays (at COMET or Mu2e), would contribute to strengthen the case for this scenario. We also illustrate how the evolution of the experimental determination of $R_{D^{(*)}}$ could be instrumental in falsifying an explanation of the anomalous $B$-meson decay data via a vector $V_1$ leptoquark.


Introduction
One of the key predictions of the Standard Model (SM) is the universality of interactions for the charged leptons of different generations. Extensive experimental observations confirm that this is indeed the case for several electroweak precision observables, as for example for Z → decays [1,2]. However, certain recent experimental measurements suggest that hints for the violation of lepton flavour universality (LFUV) might be present in a number of observables, which would thus unambiguously point towards the presence of New Physics (NP). The LFUV observables concern the flavour changing neutral current (FCNC) quark transitions b → s + − , and the charged current quark transitions b → c − ν: the former are loop-suppressed within the SM, thus providing a high sensitivity to probe NP effects; the latter can occur at the tree-level and are only subject to Cabibbo-Kobayashi-Maskawa (CKM) suppression within the SM. Among these observables, ratios of potentially LFU violating B-meson decays are of particular interest, since they are free of the theoretical hadronic uncertainties arising from the form factors, as these cancel out in the ratios. The most relevant LFUV ratios for our study are R D ( * ) (corresponding to the charged current transition b → c − ν) and R K ( * ) (corresponding to the neutral current transition b → s + − ), respectively defined as where = e, µ. A number of experimental measurements [3][4][5][6][7][8][9][10][11][12][13][14][15][16] shows deviations from the theoretical SM predictions [5,[17][18][19][20][21][22][23]. In particular, the current measurements of R D [5,11] and R D * [5,[9][10][11] respectively reveal 1.4σ and 2.5σ deviations with respect to their SM predictions [19,20,23] and, when combined, this amounts to a deviation of 3.1σ from the SM expectation [5,17,18]. In the neutral current b → s + − transitions, the measurements of R K [12,24] in the dilepton invariant mass squared bin [1.1, 6] GeV 2 show a deviation from the corresponding SM prediction [21,22] at the level of 2.5σ; for R K * , the measurements in the dilepton invariant mass squared bins q 2 ∈ [1.1, 6] GeV 2 and q 2 ∈ [0.045, 1.1] GeV 2 [13] reveal tensions with the SM expectations [21,22] with significances of 2.5σ and 2.4σ, respectively. The recent Belle collaboration results for R K * in the analogous bins [14] are consistent with both the SM and the LHCb measurements [13]. Furthermore, the LHCb measurement of R K in the bin [1.1, 6] GeV 2 has been recently updated [25], now exhibiting a 3.1 σ tension with respect to the SM prediction.
In addition to the LFUV ratios, further discrepancies with respect to the SM have also been identified in a small number of lepton flavour specific observables in b → s + − neutral current transitions -this is the case of several angular observables in both charged and neutral B 0,+ → K * µ + µ − decays (as recently reported by the LHCb collaboration [26,27]), for which tensions between observation and SM expectations lie around the 3σ level. Very recent measurements [28] of the differential branching fraction of B s → φµ + µ − decays further corroborate the picture. Moreover, LHCb recently updated [29,30] their analysis of B (s) → µ + µ − decays leading to an improved measurement of the B (s) → µ + µ − branching fractions.
In this work, our goal is to evaluate the prospects of a vector leptoquark transforming as (3, 1, 2/3) as a viable hypothesis to address the current LFUV hints, fitting in a simplfied-model approach the V 1 couplings to SM fermions, further considering how current and upcoming experimental data may strengthen or disfavour such an hypothesis. The viability of the vector leptoquark V 1 as a solution of the current LFUV hints has been explored in detail in the existing literature: in most of the existing studies, driven by an explanation of the anomalous data, only selected vector leptoquark couplings to specific quark and lepton flavours (while setting others to vanishing values) have been analysed; other studies adopt an approach driven by an ultra-violet (UV) complete model (often assisted by additional flavour symmetries) to achieve the pattern of vector leptoquark couplings preferred by the anomalous data. In contrast, we pursue a distinct avenue in this study in what concerns the underlying analysis. We rather adopt a completely data-driven approach, in which we start from a "democratic" matrix for the vector leptoquark couplings to SM quarks and leptons, without any pre-bias towards particular flavours; the existing constraints from various flavour observables, together with the anomalous LFUV data, then determine the phenomenologically allowed ranges for the vector leptoquark couplings to different flavours of quarks and leptons. We thus present a fit of the full 3 × 3 matrix of the vector leptoquark couplings taking into account all the relevant (and most stringent) measurements of various flavour observables as well as the anomalous LFUV data. As we proceed to discuss, our findings suggest that searches for a number of rare decays and transitions -conducted at Belle II and coming charged lepton flavour violation (cLFV) dedicated experiments -may help strengthening the case for V 1 models, or then contribute to exclude them as single-mediator explanations to the R K ( * ) and R D ( * ) anomalies.
Starting from a general effective theory framework, we first perform global fits taking into account the current experimental status of several observables associated with the anomalous B-meson decays, in complementary channels and kinematically interesting regions. We focus on the impact of the latest b → s data from LHCb [25][26][27][28][29][30] and how the new global fits of current flavour data favour specific classes of NP realisations.
The bulk of our work is then devoted to the study of V 1 . Numerous ultra-violet (UV) complete models for vector leptoquarks with O(TeV) mass, capable of addressing both charged and neutral current B-decay anomalies (while being consistent with the constraints from the charged lepton flavour violating decays K L → µe and K → πµe) have been proposed in the literature [82,83,87,96,102,[129][130][131][132][133][134][135]137,138], accompanied by extensive studies regarding LHC signatures and other phenomenological aspects. In this work, we pursue a distinct approach, relying on a simplified-model parametrisation of the vector leptoquark interactions with SM fermion fields, and study the impact of this scenario for a large set of observables -various leptonic and semileptonic meson decays and cLFV observables (in addition to the "anomalous" B-meson observables). The cLFV observables are very relevant not only given the current stringent constraints on the model, but also in view of the excellent (near future) projected experimental sensitivities.
Following an extensive global fit for the V 1 couplings in view of the experimental data on anomalous B-meson decay observables and data on other relevant (semileptonic) processes (meson decays and cLFV observables offering stringent constraints), we identify several tauonic modes and cLFV transitions which are likely to play a key rôle in testing the vector leptoquark scenario as a unified explanation to the B-decay anomalies. As we emphasise in this work, τ → φµ decays emerge as one of the "golden channels" for probing the V 1 hypothesis at Belle II; if the B-meson decay anomalies persist at their current level, sizeable contributions -within future experimental reach -, are also expected for B s → τ + µ − , B + → K + τ + e − , B + → K + τ + µ − and B s → φτ + µ − decays, as also noticed in [42,135]. Owing to their impressive expected future sensitivity, the upcoming cLFV experiments dedicated to searching for neutrinoless µ − e conversion in Aluminium nuclei, Mu2e and COMET, will probe a large part of the preferred V 1 parameter space via µ and e couplings to all quark generations. Furthermore, and as a direct consequence of accommodating the charged current anomalies in R D ( * ) , our study reveals that a number of b → sτ τ branching fractions are expected to be enhanced with respect to the SM (by one to two orders of magnitude), as first pointed out in [142].
Our study is organised as follows: following an EFT-based global fit in Section 2 (allowing to identify the currently best favoured NP classes of models), Section 3 is devoted to vector leptoquark realisations: in particular, and after introducing the simplified approach to this class of NP models, we perform a global fit of the V 1 couplings to SM fermions, taking into account a thorough set of flavour observables. Finally, Section 4 contains a discussion of the prospects for probing the vector leptoquark hypothesis as a solution to the anomalous B-meson decay data through several mesonic and leptonic decays to be searched for at Belle II and future cLFV-dedicated facilities. Complementary and/or detailed information relevant to the study is collected in several appendices.

Semileptonic B B B-meson decays: the impact of new LHCb data
In this section we will address the B-meson decay observables currently pointing towards a violation of LFU. Our goal is to evaluate how new recent data [25][26][27][28][29][30] has impacted the global fits carried out for an EFT approach to NP models (i.e. for generic BSM realisations), in terms of the relevant semileptonic Wilson coefficients C qq ; . We aim at investigating how well-motivated scenarios for (sets of) C qq ; -resulting from different effective operators -remain viable or have become disfavoured.
We thus begin by briefly commenting on charged current b → c ν decay data; we then proceed to discuss the experimental status of several observables associated with neutral current semileptonic B-meson decays, focusing our attention on global fits of b → s transitions, and how the status of the latter has evolved in recent months.
2.1 b → cτ ν b → cτ ν b → cτ ν data and new-physics interpretations A number of reported results from several experimental collaborations have suggested a possible violation of lepton flavour universality in the charged current decay mode B → D ( * ) ν, parametrised by the R D and R D ( * ) ratios (see Eq. (1)). The latest average values of these observables, given by the HFLAV collaboration [5], are The relevant effective Lagrangian for the charged current transitions d k → u jν − can be expressed as in which we have assumed the neutrinos to be left-handed and where, for the SM, we have For the convenient double ratios R D /R SM D and R D * /R SM D * (which combine the current experimental averages with the SM predictions), the current data can be summarised as R D /R SM D = 1.14 ± 0.10 , R D * /R SM D * = 1.14 ± 0.06, where the statistical and systematical errors have been added in quadrature.
To perform a numerical analysis of the transition B → D ( * ) τ ν (and fit the above double ratios) one further requires knowledge of the hadronic form factors which parameterise the vector, scalar and tensor current matrix elements. However, under the simplifying assumption of a non-vanishing single type of NP operator at a time -i.e. C i = 0, i ∈ {S L , S R , V L , V R , T L } -, it is possible to draw some qualitative conclusions from the approximate numerical forms for the double ratios using a heavy quark effective theory (HQET) formalism [72,[143][144][145][146][147][148].
In particular, and if one assumes that all the relevant Wilson coefficients are real, then the following qualitative observations can be readily made. The operator corresponding to C V L contains the same Lorentz structure as the SM contribution and the NP amplitude adds to the SM one, thus leading to similar enhancements to both R D and R D * , which are proportional to (1 + C V L ) 2 . In turn, this leads to similar fractional enhancements to R D /R SM D and R D * /R SM D * . Therefore, C V L is one of the most favoured choices for explaining the anomalous R D and R D * data. On the other hand, if the new physics contribution is purely a right-handed vector current (C V R type), then for a real Under such circumstances, it is then not possible to simultaneously explain both R D and R D * data. However, and as discussed in [143], this conclusion is no longer valid for a complex C V R . The scalar operators corresponding to C S L and C S R contain the pseudoscalar Dirac bilinear and therefore are not subject to helicity suppressions, leading to stringent constraints from the (relatively large) branching ratios of B c → τ ν. The tensor operator, corresponding to C T L , is subject to tensions from the recent measurement of the D * longitudinal polarisation f D * L , which is currently about 1.6 σ higher than the SM prediction and has a discriminatory power between the scalar and tensor solutions [36,40,149]. Choices based on pure right-handed operators seem to be disfavoured by LHC data [36,102]. Finally, scenarios that only present scalar contributions are in conflict with both LHC and B c → τ ν data.

Neutral current
A number of anomalies reported in b → s observables currently stand as promising hints of NP, among them those parametrised by the R K ( * ) ratios, defined in Eq. (1). The latest averages of the reported anomalous experimental data, together with the SM predictions can be expressed as [12][13][14] R LHCb where the dilepton invariant mass squared bin (in GeV 2 ) is identified by the associated subscripts. Further anomalies have also been reported in the neutral current decay modes of B-mesons for semileptonic final states including muon pairs 1 . Among them, one concerns the observable Φ ≡ dBR(B s → φµµ)/dq 2 in the bin q 2 ∈ [1, 6] GeV 2 [15], presently exhibiting a tension with the SM prediction around 3σ. Further discrepancies with respect to the SM, typically at the 3σ level, have also emerged in relation to the angular observables. In particular, this is the case of P 5 in B → K * + − processes: the results from the LHCb collaboration for P 5 regarding muon final states (B → K * µ + µ − decays) reveal a discrepancy with respect to the SM [150,151]; the Belle collaboration [16,152] reported that P 5 results for electrons show a better agreement with theoretical SM expectations than those for muons. More recently, similar measurements have also been reported by the ATLAS [153] and CMS [154] collaborations. The 2015 LHCb results [151] and the ATLAS result [153] for P 5 in the low dimuon invariant mass-squared range, q 2 ∈ [4, 6] GeV 2 , indicate a ≈ 3.3σ discrepancy with respect to the SM prediction [155]. Belle results corroborate the latter findings, showing a deviation of 2.6σ from the SM expectation in the bin q 2 ∈ [4, 8] GeV 2 [16]. The reported CMS measurement (possibly as a consequence of insufficient statistics) is still consistent with the SM expectation within 1σ [154]. Among the angular observables it is important to stress that F L , P 4 , P 5 and P 8 have been a driving force in the evolution of the global fits. Very recently, the LHCb collaboration has updated the results for the angular observables relying on 4.7 fb −1 of data [26,27]: local discrepancies of 2.5σ and 2.9σ, respectively in the bins q 2 ∈ [4, 6] GeV 2 and q 2 ∈ [6, 8] GeV 2 GeV 2 , were reported. While these lepton flavour dependent observables are also sensitive to the presence of NP [156][157][158][159][160], they are nevertheless subject to hadronic uncertainties (for example form factors, power corrections and charm resonances [161][162][163][164][165][166][167][168][169][170][171][172]) contrary to the LFUV ratios, which are in general free of the latter sources of uncertainty.
A way to consistently analyse the aforementioned anomalous experimental data is to adopt the "effective approach", in which all possible short-distance NP effects are encoded in the Wilson coefficients related to a complete EFT basis. Within a weak effective theory (WET), the effective Lagrangian for a general d j → d i − + transition can be expressed as [173][174][175][176][177][178] with V ij denoting the CKM matrix and in which the relevant operators are defined as where the primed operators O 7,9,10,S,P correspond to the exchange P L ↔ P R . Given the above WET parametrisation, the first question to address concerns the set(s) of Wilson coefficients seemingly preferred by the anomalous experimental data, which then leads to the identification of possible phenomenological candidates, and ultimately to the construction of UV complete extensions of the SM.
Let us then first proceed to obtain model-independent fits for different possible new physics scenarios (in terms of non-vanishing contributions to one or several Wilson coefficients in the FCNC b → s transitions), with a particular emphasis on the impact of the recent data from the LHCb collaboration [25][26][27][28][29][30].
We thus begin by performing a fit of the data on angular distributions, differential branching fractions and the LFUV ratios R K ( * ) -excluding the new measurements of LHCb [25][26][27][28][29][30] -to establish a baseline to study the impact of the new measurements from LHCb. The underlying methodology for the fit as well as the details of the statistical methods are described in Appendix A; the specific bins of observables and datasets used for the fit are presented in Appendix B.1.
Using these data sets, one can already infer a qualitative behaviour of the fits in terms of the Wilson coefficients (allowing to identify favoured NP "scenarios"). While NP contributions exclusively to C bsµµ 9 already give a very good fit when compared to the SM [22, 31-35, 38, 45, 121, 179-183], most realistic NP models considered to explain the tensions also generate non-zero contributions to other Wilson coefficients, either by construction or then through operator mixings which occur when renormalisation group (RG) running effects (from the NP mass scale to the observable scale) are taken into account. In particular, the SU (2) L conserving scenario ∆C bsµµ 9 = −∆C bsµµ 10 has received a considerable attention in recent years, as it provides a very good fit to the data. However, following the improvement in the measurement of R K in 2019 [12], with relatively smaller experimental uncertainties, the preference has been slightly shifted to more involved scenarios, calling upon a larger number of non-vanishing Wilson coefficients. This becomes manifest through tensions between the individual fits for the LFUV ratios, and for the lepton flavour dependent observables (Φ ≡ dBR(B s → φµµ)/dq 2 and the angular observable P 5 in B → K * + − ), under the hypotheses of ∆C bsµµ 9 = −∆C bsµµ 10 (and C bsµµ 9 ) NP "scenario(s)". Therefore, if the anomalous data for the lepton flavour dependent observables is not due to statistical fluctuations or to long-distance effects, then such tensions suggest non-trivial NP contributions in other lepton flavours, or in distinct Wilson coefficients. For example, in Ref. [34], it was reported that LFUV contributions in b → se + e − , in addition to a minimal "scenario" (i.e. ) can ease such tensions, further improving the overall fit to data with respect to the SM. On the other hand, in [184] it was observed that if one considers a LFUV scenario which only affects muons in conjunction with a non-vanishing LFU NP contribution (i.e. with equal contributions to e, µ, and τ ), then the anomalous LFUV ratio data can be explained by the LFUV in the muon sector (with sub-leading interferences with LFU NP contributions), while the lepton flavour dependent observables can be fitted combining LFUV and LFU NP, with improved agreement with respect to the overall data. Therefore, in our global fit scenarios, we include the above two interesting possibilities, comparing individual and combined effects.
As can be seen in Table 1, we find that LFU contributions to C bsµµ 9 and C bsee 9 (∆C univ.

9
, corresponding to the last two blocks in the table) are able to significantly improve the model-independent fits. It is interesting to note that these scenarios naturally arise in many simple models attempting a combined explanation of b → s data together with the anomalous charged current data on b → cτ ν. In particular, a sizeable contribution to the charged current Wilson coefficients to explain b → cτ ν calls upon large τ -couplings, which in turn generate a sizeable C bsτ τ 9 . Through RG operator mixing effects [185] (evolution from NP scale to the observable scale), the C bsτ τ 9 contribution leads to a LFU contribution for both C bsµµ 9 and C bsee 9 . We notice that this LFUV and LFU NP combined "scenario" is of relevance for SU (2) L -singlet vector-leptoquark models since, due to the SU (2) L representation, the charged current couplings in b → cτ ν transitions are identical to the ones appearing in the neutral current b → s transitions (up to CKM elements in the effective Wilson coefficients).   Table 1: Well-motivated NP "scenarios" (sets of Wilson coefficients) and corresponding fits to the data on angular distributions, differential branching fractions and the LFUV ratios R K ( * ) (not including the new measurements of LHCb [25][26][27][28][29][30]). The SM p-value is found to be ∼ 6.8%.

Fits of
To estimate the impact of the recent measurements of LHCb [25][26][27][28][29][30], we repeat the Wilson coefficient fit of the previous subsection, but now taking into account the new LHCb data [25][26][27][28][29][30]. We further take into account the recently improved limits on BR(B 0 → e + e − ) and BR(B s → e + e − ), [186], and a recently improved measurement of the angular observables in B 0 → K * e + e − at very low q 2 , as reported by the LHCb collaboration [187].
In the fits carried out for the older data sets, the scenario ∆C bsµµ 9 = −∆C bsµµ 10 had a larger p-value than the fit which only included NP contributions to C bsµµ 9 ; as can be seen from Table 2, the situation is now reversed upon inclusion of the new LHCb data. Furthermore, ∆C bsµµ 9 = −∆C bsµµ 10 arguably provided an equally good fit (c.f. Table 1) to the data compared to the hypotheses which included a universal contribution to C bs 9 in addition to the (V − A) contribution, whereas now the hypotheses with a universal contribution are clearly preferred. In Fig. 1 we present the likelihood contours for the "pre-2020" and recent data, around the corresponding best-fit points, where it can be seen that a non-vanishing universal contribution to C 9 is now preferred at around ∼ 3σ. Although the position of the best fit point is only slightly changed (from the former diamond to the current star), the new measurement leads to an improved precision for the model-independent fits. This is manifest from the comparison of the likelihood contours belonging to either dataset (regions delimited by dashed or solid lines, respectively in association with "pre-2020" and full data).
This renders models that attempt a combined explanation of the charged and neutral current Bdecay anomalies, especially single-mediator scenarios, even more preferable. Following this section's discussion, in the remainder of our study we will focus on such NP realisations, in particular extensions of the SM via single left-handed vector fields (in our case, a V 1 vector leptoquark), which provide the best fit among all the possibilities for single-mediator NP scenarios [36,40,149].   , corresponding to the scenario with the largest p-value (see Tables 1 and 2). The shaded regions (delimited by full lines) correspond to the 1, 2 (3)σ regions around the best-fit point including the recent data; the dashed lines denote the same likelihood contours, without the inclusion of the recent LHCb measurement. In addition to the angular observables and R K ( * ) , the "global" contour (green regions) includes all other b → s data as listed in Appendix B.1. The black pentagon denotes the SM-value, while the star (diamond) denotes the best-fit point to the current (old) data.

Implications for V 1 vector leptoquark solutions
Among the many possible SM extensions including leptoquarks, in what follows we focus on vector leptoquark (V 1 ) scenarios. This possibility has received increasing attention in the literature, as it is currently the only single-leptoquark construction that successfully offers a simultaneous solution to both charged and neutral current B-meson decay anomalies [82,87,96,102,[129][130][131][132][133][134][135][136][137][138]. As highlighted following the updated global fits carried out in the previous section, the vector leptoquark hypothesis belongs to the class of NP "scenarios" most favoured by current data.
However, and in order to account for experimental data, V 1 should have non-universal couplings to quarks and leptons, and the latter can be realised in a number of ways. The most minimal possible scenario relies in the assumption that the vector leptoquark is an elementary gauge boson 2 , associated to a non-abelian gauge group extension of the SM, under which the SM fermion generations are universally charged; in the unbroken phase of the underlying extended gauge group, the leptoquark gauge couplings also remain universal. Despite its simplicity, this scenario is challenged by constraints from the cLFV decays K L → µe and K → πµe: current limits force the mass of such a vector leptoquark to be very heavy, m V ≥ 100 TeV for O(1) couplings [190][191][192][193][194][195], and thus excessively heavy to account for both the charged and neutral current B-meson decay anomalies. In order to understand this, notice that while V 1 has a universal coupling to SM fermions in the unbroken phase, after SU (2) Lbreaking a potential misalignment of the quark and lepton eigenstates is generated, leading to LFUviolating V 1 couplings. Given the constraints from τ decays, the cν coupling generated from bτ through CKM mixing is not sufficiently large to account for R D ( * ) data [196]. On the other hand, for a maximal cν coupling (with the neutrino flavour in cν different from ν τ ) generated by d i µ and d i e couplings, important constraints arise from R K ( * ) data for i = 2, 3, and from kaon decays for i = 1. Moreover, the cν coupling induced by dτ is heavily CKM-suppressed. Therefore, the only viable possibility is to maximise both bτ and sτ couplings, which in turn will induce large couplings between the first two generations of quarks and leptons (given the unitarity of the post-SU (2) L -breaking mixing matrix), thus implying excessive contributions to cLFV.
A possible way to circumvent the above mentioned constraints is to introduce three "generations" of vector leptoquarks, belonging to an identical number of copies of the extended gauge group (e.g. Pati-Salam model based on the gauge group [SU (4) , with each copy acting on a single SM fermion generation (subject to mixing with additional vector-like fermions), with the largest leptoquark-fermion couplings in association with the third family [131]. Another possibility to lower the vector leptoquark mass relies in an extended gauge group, [83]. This leads to an approximate U (2) flavour symmetry, which is softly broken by new vector-like fermions, thus allowing to obtain the desired non-universality in the leptoquark couplings. An alternative simplified-model framework, without the need to specify an explicit extended gauge group, was pursued in [96]: working under a single vector leptoquark hypothesis, an effective non-unitary mixing between SM leptons and new vector-like leptons was used to account for the LFUV structure required to simultaneously explain both the charged and the neutral current B-meson decay anomalies.
Irrespective of the actual NP model including (not excessively heavy) vector leptoquarks, the effects can be understood in terms of contributions to the Wilson coefficients. Following the discussion of the previous section (see Tables 1 and 2), in order to achieve the preferred contributions for the Wilson coefficients, C bsµµ 9 = −C bsµµ 10 and a universal ∆C univ. 9 , scenarios in which V 1 couples at the tree level through a left handed (V − A) current to muons (as well as to down-type quark flavours b and s) appear to be favoured the most by the global fits. A nonvanishing ∆C bsee 9 = −∆C bsee 10 along with C bsµµ 9 = −C bsµµ 10 and a universal ∆C univ. 9 also provides a reasonable fit but such hypotheses are subject to stringent constraints from cLFV processes. Furthermore, and in order to also address the charged current data (R D ( * ) ), sizeable tree-level τ couplings to second and third generation quarks must also be present, and these induce new contributions to the C V L Wilson coefficient. Such large V 1 − τ couplings to second and third generation quarks further lead to a large C bsτ τ 9(10) which then feeds into the muon and electron counterparts (in a universal way) through RG running 3 . A simplified-model parametrisation of the vector leptoquark couplings allows not only to perform global fits, but also to understand the phenomenological implications of the relevant flavour structure, which is paramount to establish the current viability of the model, and its prospects for future testability. In this section, we thus pursue this approach not only regarding the "anomalous" B-meson observables, but also in what concerns the impact of this BSM construction for a large set of observables (various flavour violating meson decays and cLFV modes) -relevant in terms of constraints on the model, or then offering excellent prospects of observation in the near future.

A simplified-model parametrisation of vector leptoquark V 1 couplings
As mentioned before, in our study we will focus on SM extensions via a vector leptoquark V 1 , arising from an unspecified gauge extension of the SM. The new vector transforms as ( For simplicity, and due to the absence of hints in the data suggesting the presence of right-handed couplings 4 , we will exclusively focus on left-handed leptoquark currents. In the mass basis we consider a simplified-model Lagrangian concerning the effective coupling of V 1 with the SM fermions, given by in which K ij L are effective couplings which are in general complex and non-universal, V denotes the CKM matrix and U P ≡ U † L U ν L is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) leptonic mixing matrix. We notice that the above parametrisation is valid irrespective of the underlying mechanism responsible for the generation of the effective nonuniversality in the vector leptoquark couplings (see, e.g. [83,87,96,102,131,[136][137][138]). For the sake of simplicity we will further assume that the couplings V 1 − − q are real. The conclusions drawn here should thus hold for generic constructions with real effective couplings and negligible right-handed currents (consistent with zero at the ∼ 1σ level from the EFT fits to the b → s data [31]; we recall that this corresponds to negligible C 7, 9, 10 Wilson coefficients).
For the general vector leptoquark scenario under consideration, the most relevant tree-level Wilson coefficients for b → s transitions and R D ( * ) observable are given by [197] C ij; Variants of the above coefficients (depending on the flavour indices) are responsible for the leading contributions to most of the b → s and R D ( * ) observables relevant for the fit. In addition, there are several other observables such as leptonic and semileptonic meson decays, as well as cLFV leptonic decays, which are important for the analysis. The expressions for the branching fractions can be found in Appendix C. Before proceeding to the description of the global fit, some remarks are in order concerning the evaluation of the latter observables. We first notice that the vector leptoquark coupling parameters are matched with the Wilson coefficients at the leptoquark mass scale, m V ; the latter are subsequently run down to the b-quark mass scale, or to the scale of any other process (observable) considered in the analysis. Therefore, and even though some of the relevant Wilson coefficients are vanishing at the scale of the matching of the EFT to our simplified effective leptoquark model (i.e. at m V ), they can be generated from operator mixing during RG running to the scale of a given observable. In particular, in our fits we take into account all running effects using the wilson package [198] in association with the flavio package [199]. We recall that the non-trivial effective V 1 couplings can potentially induce new contributions to cLFV observables such as radiative decays i → j γ and 3-body decays i → 3 j at loop level, and neutrinoless µ − e conversion in nuclei (at tree-level). In view of the very good current experimental sensitivity, these observables will provide some of the most stringent constraints on the V 1 couplings to SM fermions; as already mentioned, the expected improvements on the future sensitivities offer the possibility to further probe the vector leptoquark couplings. Another important point worth noting is that although the cLFV radiative decays occur at loop level, the associated anapole contributions to vector operators can lead to sizable contributions to neutrinoless µ − e conversion and µ → 3e, with a magnitude comparable to the tree level contributions or, in some cases, even accounting for the dominant contribution. In addition, we find that the dipole operators also significantly contribute to radiative decays and to neutrinoless µ − e conversion.
We emphasise that the one-loop dipole and anapole contributions from the exchange of a vector boson generically diverge, and a UV complete framework (with a consistent gauge symmetry breaking pattern) is thus necessary to obtain a finite result in a gauge independent way. Therefore, to reliably evaluate such observables in the context of vector leptoquark exchanges we have chosen to work in the Feynman gauge, including the necessary contributions from the Goldstone modes. We have thus made the minimal working assumption that the vector leptoquark originates from the breaking of a gauge extension of the SM, which gives rise to a would-be Goldstone boson degree of freedom, subsequently absorbed by the massive vector leptoquark. We include this Goldstone mode (degenerate in mass with V 1 ) to obtain the gauge invariant (finite) form factors for the relevant dipole and anapole contributions. Furthermore, to keep our results as general as possible, we do not include any effects due to an extended scalar sector (possibly necessary to implement the breaking of the extended gauge symmetry) nor from new gauge bosons (which might arise due to the breaking of an extended gauge symmetry); we work under the assumption that, should these states be present, they only have flavour conserving couplings 5 .
Finally, and should the model encompass additional vector-like fermions, as is the case for several vector leptoquark realisations [83,96,131], one must also consider the impact of such new states for electroweak precision observables, as for instance the constraints on the Z boson LFU ratios and cLFV decay modes, as emphasised in [96].

Towards a global fit of the vector leptoquark V 1 flavour structure
We are now ready to carry out a comprehensive fit of the relevant couplings of the vector leptoquark to the different generations of SM fermions. Relying on the above simplified-model parametrisation, our goal is thus to constrain the entries of the matrix K L (see Eq. (7)). Under the assumption that the relevant couplings are real, a total of nine free parameters will thus be subject to a large number of constraints stemming from data on several SM-allowed leptonic and semileptonic meson decays, SMforbidden cLFV transitions and decays, as well as from an explanation of the (anomalous) observables in the b → s and b → cτ ν systems.
Other than studying the contributions of the vector leptoquark in the "anomalous" channels, we aim to estimate the favoured ranges of all of its couplings to SM fermions. Consequently, we include a large number of additional observables into the likelihoods. Since most processes only constrain a product of at least two distinct leptoquark couplings, a successful strategy is to include an extensive set of processes, thus allowing to constrain distinct combinations of couplings (as many as possible).
In addition to the b → c ν transitions, we also include certain b → u ν decays such as B 0 → πτ ν, B + → τ ν and B + → µν, which are listed in Table B.5. In many leptoquark models B → K ( * ) νν decays provide very stringent constraints. However this is not the case for V 1 vector leptoquarks, due to the SU (2) L -structure: the relevant operators for B → K ( * ) νν transitions are absent at the tree-level, and are only induced at higher order, thus leading to weaker constraints. Due to the leading operator being generated at the loop level, a non-linear combination of leptoquark couplings is constrained by this process. Thus, despite the loop suppression, we include B → Kνν in the likelihoods, and use the data obtained by Belle [218,219] and BaBar [220,221].
To constrain combinations of first and second generation couplings, we further include a large number of binned and unbinned leptonic and semileptonic charged current D meson decays, charged and neutral current kaon decays and SM allowed τ -lepton decays. The observables and corresponding data-sets can be found in Appendix B.3 and are listed in Tables B.6 through B.8.
Finally, cLFV processes impose severe constraints on the parameter space of vector leptoquark couplings; in particular neutrinoless µ − e conversion in nuclei and the decay K L → e ± µ ∓ provide some of the most stringent constraints for vector leptoquark couplings to the first two generations of leptons [96]. In Table 3 we present the current experimental bounds and future sensitivities for various cLFV observables yielding relevant constraints to our analysis. Depending on the fit set-up, either only a few, or then all of these observables are included in the global likelihood, as explicitly mentioned in the following paragraphs.   Results for the simplified-model fit of the V 1 couplings Firstly, it is important to emphasise that in our analysis we consider all the entries in the K L coupling matrix as (real) free parameters to be determined by the fit. For the leptoquark mass we choose three benchmark-points, m V 1 ∈ [1.5, 2.5, 3.5] TeV, which allow to illustrate most of the vector leptoquark mass range of interest, while respecting the current bounds from direct searches at colliders [237][238][239][240][241][242][243][244][245]. In particular, notice that masses significantly heavier than a few TeVs preclude a successful explanation of the charged current anomalies, R D ( * ) . For each mass benchmark point we thus obtain best-fit points corresponding to a SM pull around ∼ 6.4 σ (with respect to the global likelihood including all lepton flavour conserving observables). In Fig. 2, we present the results of a random scan around the best-fit points for the vector leptoquark scenario here considered, in the plane spanned by two of the most constraining cLFV observables, CR(µ − e, N) and BR(K L → e ± µ ∓ ). The sample points are drawn from the posterior frequency distributions of the leptoquark couplings, following Markov Chain Monte Carlo (MCMC) simulations, as described in Appendix A. It can be easily seen that for the three mass benchmark choices (corresponding to the different colours in the plot) most of the randomly sampled points are excluded by the strong cLFV constraints. Although the involved couplings are compatible with 0, the constraints on first generation couplings derived from lepton flavour conserving low-energy data (as listed in Appendix B) are considerably weaker than those from LFV processes. This leads to several "flat directions" in the likelihood. The strongest LFV constraints are from CR(µ − e, Au) and BR(K L → e ± µ ∓ ), while other LFV constraints on second and third generation couplings are weaker, or on par with constraints from lepton flavour conserving low-energy data. Therefore, we redefine the strategy of the global fit, and now directly include the upper bounds from CR(µ − e, Au) and BR(K L → e ± µ ∓ ) as inputs in the fitting procedure for the vector leptoquark couplings.
The inclusion of the current upper limits on the observables CR(µ−e, Au) and BR(K L → e ± µ ∓ ) as input to the fit will consequently shift the best-fit point towards a lower cLFV prediction, also leading to a slightly lower SM pull. However, we find this to be a good compromise in order to identify regimes in the parameter space not yet disfavoured by the current cLFV data. In fact, and since CR(µ − e, Au) and BR(K L → e ± µ ∓ ) are indeed two of the most constraining cLFV observables, once the bounds on the latter observables are respected, most of the sample points will be naturally in agreement with current bounds on most of other cLFV observables (this is a consequence of correlations with other cLFV µ − e transitions; processes involving τ -leptons are comparatively less constraining).
In Table 4 we present our results for the new fits with their corresponding SM pulls. As can be verified, the SM pull is lower, reduced from ∼ 6.4σ to ∼ 5.8σ, of which the contributions to the total χ 2 stemming from the charged current b → c ν transitions amounts to ∼ 1.5σ, whereas the contributions from the neutral current b → s transitions amounts to ∼ 4.3σ. Furthermore, we show tentative 90% ranges of the posterior (coupling) distributions, obtained by sampling the global b → s likelihood using MCMC. The ranges, derived from the histograms of the posterior distributions, are taken as symmetric intervals between the 5 th and 95 th percentiles (cf. Appendix A). We notice here that the vector leptoquark coupling to the first generation SM fermions are consistent with zero, which is an assumption often invoked in literature for simplified analyses. For second and third generation couplings, the quoted ranges of the corresponding fits are in fair agreement with the (order of magnitude) results for the benchmark ranges of second-and third generation couplings quoted in the literature, e.g. [135,246]. However, given the differences in the coupling parametrisation choices and underlying statistical treatment, the results are not directly comparable.  Table 4: Results of the fits including the current experimental bounds on CR(µ−e, Au) and BR(K L → e ± µ ∓ ) in the likelihood: best fit points and symmetric 90% ranges (see Appendix A for details) of K ij L . The SM pull is reduced from ∼ 6.4σ to ∼ 5.8σ.
Upon inclusion of the current cLFV constraints we find that the shape of the global likelihood consequently enforces small vector leptoquark couplings to the first two generations of charged leptons, leading to predictions consistent with experimental data. This thus allows to sample the global likelihood (in terms of the leptoquark couplings) via MCMC techniques (as described in Appendix A). The posterior distributions of the leptoquark couplings are then used to compute predictions for B-meson decays into final states containing τ -leptons, and several cLFV observables (including tau decays). This is presented in Fig. 3 where, for each observable, we depict the current experimental bounds and future sensitivities, the SM predictions (when relevant), as well as the predictions for the three vector leptoquark mass benchmark points -corresponding to the vertical coloured lines. The dashed lines describe predictions of observables involving only couplings compatible with vanishing values and thus their top edge corresponds to a 90% upper limit, while no lower limit should be implied.
As can be seen from Fig. 3, a large part of the currently allowed parameter space in the eµ channel (for the three leptoquark mass benchmark points) will be probed by the upcoming experiments dedicated to searching for neutrinoless µ−e conversion in Aluminium nuclei, Mu2e and COMET, owing to the expected increase in sensitivity. In the case of future non-observation of this process, this will lead to strongly improved constraints on the V 1 couplings to first first generation fermions.
Moreover, the sensitivity of the lepton flavour violating process τ → φµ is expected to be improved by over an order of magnitude at the Belle II experiment, which will allow probing a large region of the parameter space associated with the µτ channel. A priori, and as can be seen from Fig. 3, under the current vector leptoquark hypothesis, τ → φµ decays have very strong prospects of being observed at Belle II. Conversely, should such a mode not be observed at Belle II, then the s − µ and s − τ couplings of the vector leptoquark will be tightly constrained. As a consequence, it might prove extremely challenging to simultaneously address the anomalous neutral and charged current data within the current model.

Impact of future experiments: Belle II and cLFV searches
Following the overview of the vector leptoquark couplings conducted in the previous section, we now proceed to investigate how our working hypothesis can be effectively probed by the coming future experiments, especially Belle II and cLFV-dedicated facilities.  Assuming that the above experiments return only negative search results for the most promising modes, we then evaluate how the current V 1 hypothesis would still stand as a viable explanation for the LFUV B-meson decay anomalies.

Probing the vector leptoquark V 1 at coming experiments
Concerning the quest for LFUV in b → s + − decays, Belle II is expected to achieve a very high sensitivity for both muon and electron modes, leading to very precise measurements for the ratios R K and R K * , with the potential to confirm the anomalous LHCb data (if the latter is due to NP effects) [225]. In what concerns B-meson decays to τ + τ − final states, Belle II will also provide the first in-depth experimental exploration of these modes. Notice that the latter remain a comparatively less explored set of observables, with relatively weak bounds on the few modes already being searched for: for example, current bounds on BR(B 0 → τ + τ − ) < 1.3 × 10 −3 from LHCb [247] and BR(B s → τ + τ − ) < 2.25 × 10 −3 from Babar [248] are orders of magnitude weaker than the SM predictions. For the purely leptonic decays, the most recent SM computations now include next-to-leading order (NLO) electroweak corrections and next-to-NLO QCD corrections [249][250][251], Within the SM, the exclusive semileptonic decays of B-mesons to τ + τ − final states have been studied by several groups: the modes B → K * τ + τ − and B s → φτ + τ − have been computed 6 in [254][255][256].
As already discussed in Section 3, sizeable b − τ and s − τ couplings are necessary to explain the charged current anomalous data on R D ( * ) ; if R D ( * ) anomalies are indeed due to NP then one expects a significant enhancement of the rates of b → sτ + τ − processes, up to three orders of magnitude from the SM predictions [70,78,130,142]. This renders searches for b → sτ + τ − modes extremely interesting probes of vector leptoquark models aiming at explaining anomalous LFUV data. Although the LHCb programme includes searches for B → K ( * ) τ + τ − and B s → φτ + τ − modes, being an e + e − experiment Belle II is expected to be more efficient than the LHCb in reconstructing B to tau-lepton decays, since many of these modes require reconstructing additional tracks originating from the final state mesons (K, K * or φ). Therefore, b → sτ + τ − observables will be among the "golden modes" aiming at probing the vector leptoquark hypothesis at Belle II.
In Fig. 4 we present the predictions for several leptonic and semileptonic B (s) to τ + τ − decays, as arising in the present vector leptoquark scenario. We display the results for three benchmark leptoquark masses (coloured vertical bars, corresponding to m V =1.5 TeV, 2.5 TeV and 3.5 TeV), together with the current limits and the future projected sensitivity from Belle II, and the corresponding SM predictions.
As can be clearly observed from Fig. 4, all b → sτ τ branching fractions are enhanced with respect to the SM (typically by one to two orders of magnitude). This is a direct consequence of accommodating the charged current anomalies (i.e. R D ( * ) ), as these call upon sizeable b − τ and s(c) − τ couplings. The decay B 0 → τ + τ − is subject to a milder enhancement due to having the d − τ coupling already constrained by other observables.
Tau-lepton decays offer powerful probes of vector leptoquark models. The Belle experiment has searched for 46 distinct cLFV τ decay modes, using almost its entire data sample of approximately 1000 fb −1 ; no evidence for cLFV decays was found, but new 90% C.L. upper limits on the branching fractions were set, at a level of around O(10 −8 ). At Belle II, if on the one hand the higher beam-induced background will render these searches more challenging, on the other hand its impressive luminosity will allow to significantly ameliorate the sensitivities to these modes. As much as 45 billion τ pairs (in the full dataset) are expected to be produced in e + e − collisions at Belle II, clearly providing very bright prospects for cLFV tau decay searches. The Belle II experiment is thus expected to improve dBR dq 2 [15,22] ( the sensitivities of the various cLFV decays by more than one order of magnitude, reaching a level of O(10 −9 − 10 −10 ). In Fig. 5 we present the predictions of the vector leptoquark scenario for various cLFV tau decay modes which are programmed to be searched for at the Belle II experiment. It is interesting to note that among the various observables, the τ → φµ decay emerges as the most promising one to probe the vector leptoquark hypothesis -another "golden mode".
The Belle II experiment will also search for a number of cLFV leptonic and semileptonic Bmeson decays (some into final state τ s). In Fig. 6 we present our findings for these cLFV processes. In the context of the present vector leptoquark model, one thus expects sizeable contributions for B s → τ + µ − , B + → K + τ + e − , B + → K + τ + µ − and B s → φτ + µ − (for the different benchmark masses considered), close to current bounds, and clearly within reach of future sensitivities 7 . Together with the decay channels identified following the results displayed in Fig. 4, these cLFV modes appear particularly promising to observe a signal of a vector leptoquark NP scenario explaining the B-meson decay anomalies.

Impact of future negative searches
A final point to be addressed concerns the impact of future null results from Belle II and other experiments searching for cLFV: if no cLFV signal is found, and no enhancement of B-meson decay rates is observed, to which extent will this affect the prospects of a vector leptoquark hypothesis as a viable explanation of the B-meson decay anomalies? To assess the implication of such a scenario we re-conduct the fit whose results were summarised in Table 4, now including the projected future sensitivities from Belle II and cLFV-dedicated experiments (COMET, Mu2e, MEG II and Mu3e). 7 Notice that the rates for Bs decays into φτ − µ + are typically less enhanced than those for the (opposite charge) φτ + µ − mode: this is a consequence of the leptoquark couplings involved, with the combination K 22 L K 33 L (entering the former) in general smaller than K 23 L K 32 L (appearing in the latter), as can be inferred, for example, from Table 4. The 90% ranges are obtained from sampling points at the around the best-fit point. Line and colour coding as in Fig. 3. Recall that the Belle II observables taken into account in this fit are listed in Appendix B (Table B.9), with the future sensitivities always corresponding to the assumption of the full anticipated luminosity of 50 ab −1 ; the future sensitivities for the cLFV dedicated experiments have been summarised in the first part of Table 3.
The results of this new fit (corresponding to null results in the several "golden modes" previously discussed) are presented in Table 5. A comparison of these results with those of Table 4 suggests that all leptoquark couplings would be well constrained (with the exception of the d − τ one). We again notice here that the vector leptoquark coupling to the first generation SM fermions remain consistent with zero.  One can now re-project the new fit results onto the plane of the anomalous B-meson decay observables, by randomly sampling around the best fit points presented in Table 5. For the V 1 scenario under consideration, the strongest impact of a non-observation of cLFV processes and non-enhanced rates for B-meson decays to τ + τ − final states occurs for the fit of the charged current anomalies R D and R D * . This is a consequence of having significantly stronger constraints on the vector leptoquark couplings to τ -leptons following the negative search results from Belle II and future cLFV experiments, and will render V 1 less efficient in contributing to both R D ( * ) .

TeV
We present in Fig. 7 the different likelihood contours and leptoquark predictions, for different benchmark masses 8 and fit set-ups, as well as best-fit points for the distinct experimental scenarios. The impact for the b → c ν fit can be observed in the R D − R D * plane depicted in Fig. 7, as the preferred "region" (orange cross) is pulled towards the SM prediction, and away from the current experimental best fit point (red circle).
Notice however that potential negative results from Belle II and future cLFV experiments do not significantly affect the fit to anomalous b → s observables.
The above discussion clearly emphasises the key rôle played by Belle II and future cLFV experiments in probing the vector leptoquark scenario as a unified explanation to the B-decay anomalies, especially in view of a new determination of R D ( * ) (central value and associated uncertainties). Scenarios can be envisaged in which future experimental data corroborates current R D ( * ) values (no change in the central value, corresponding to the red "dot" in Fig. 7), but accompanied by a reduction of the associated errors (implying tighter likelihood contours): this could then potentially contribute to disfavour V 1 as a viable explanation to the charged current B-meson decay anomalies. However, if future Belle II data (dashed contours in Fig. 7) evolves along current Belle data, vector leptoquarks would still remain exceptional candidates to explain the B-meson decay anomalies, while avoiding detection in cLFV processes in the future.

Concluding remarks
Being a well-motivated new physics candidate, leptoquark extensions of the SM have been increasingly investigated, in view of their potential for a simple, minimalistic scenario to explain the current hints of LFUV arising from B-meson decay data. Vector leptoquarks transforming as ( Figure 7: Likelihood contours and vector leptoquark predictions for R D and R D * . Red regions correspond to different likelihood contours obtained from a naïve combination of the experimental likelihoods. The blue cross denotes the predictions at the best-fit point to current LFV data. The orange cross denotes the predictions at the best-fit point with assumed null results of LFV processes at Belle II, Mu2e and COMET. The black cross denotes the SM prediction [5]. The green dashed contour line describes the naïve extrapolation of the current combination of Belle data [6,10,11] to the anticipated future precision of the Belle II experiment, while the purple dashed contour line is a naïve combination of the Belle II projection with the current data. particularly appealing, as they offer a simultaneous explanation for both charged and neutral current B-meson decay anomalies, parametrised by the R K ( * ) and R D ( * ) observables. In our work, we have thus investigated how minimal constructions, containing the vector leptoquark V 1 , successfully account for the anomalies in both R K ( * ) and R D ( * ) . Leading to our study, and relying on an EFT approach, we first presented results of global fits, which allowed to assess the impact of the most recent LHCb data in identifying the most favoured generic classes of NP realisations (in terms of new contributions to the relevant Wilson coefficients). Our findings suggest that scenarios in which a universal contribution to C bs 9 is present -in addition to the (V − A) contribution -become increasingly preferred (at a ∼ 3σ level).
In our study we have taken a phenomenological approach for the couplings of the vector leptoquark to SM quarks and leptons: we emphasise that starting from a completely general simplified-model parametrisation, we presented a fit for the full 3 × 3 matrix (K ij L ) encoding the V 1 q couplings, taking into account various relevant flavour observables and the anomalous LFUV data. In addition to providing a better guidance towards possible UV completions of vector leptoquark scenarios capable of addressing the LFUV anomalies, this approach can also reveal interesting prospects for observables which can be potentially used to probe the underlying vector leptoquark hypothesis (we notice that many of the latter observables can be missed in analyses with a priori vanishing couplings to the first generation of SM fermions. Relying on this alternative formalism for the phenomenological fitting of the vector leptoquark couplings, we thoroughly investigated the impact of such a NP scenario: we considered the prospects for an extensive array of observables, including (in addition to the anomalous B-meson decay observables) leptonic cLFV transitions, several B decay modes to final states including τ + τ − pairs, flavour violating τ decays as well as cLFV (semi)leptonic decays of B-mesons. In view of the excellent experimental prospects, we have investigated several very promising "golden modes" to (indirectly) test the V 1 scenario. Among these channels one finds τ → φµ decays, b → sτ τ and b → sτ µ transitions, as well as µ − e conversion in nuclei. These modes, searched for at Belle II and coming cLFV experiments, will play a crucial rôle in testing the vector leptoquark hypothesis as a single explanation to the R K ( * ) and R D ( * ) anomalies.
As we have discussed, the confirmation of LFUV in B-meson decays, (strongly) enhanced rates for B-meson decays to τ + τ − final states, as well as an observation of cLFV transitions in certain channels (by itself a massive discovery!), would all contribute to substantiate a vector leptoquark NP scenarioalthough some of the latter signals could indeed arise from other BSM constructions. Conversely, the non-observation of such signals at Belle II and future cLFV experiments has the potential to falsify the vector leptoquark scenario as a solution to the anomalous R D ( * ) data, if the latter anomaly persists in future measurements with reduced uncertainty (without significant changes in the central values). Should this be the case, and although NP models containing vector leptoquarks could still address the neutral current B-decay anomalies (i.e. R K ( * ) ), a common explanation of both sets of anomalies would be certainly more challenging.
The coming years clearly offer rich and promising experimental prospects to test one of simplest -yet successful -new physics constructions that allows explaining both the LFUV B-meson decay anomalies.

A Statistical treatment and fits
A proper statistical treatment of the experimental data and of the theoretical uncertainties is imperative for a precision analysis of flavour observables. In general, the goal is to find a set of theoretical predictions for the observables of interest (O th i ) which agrees best with the experimental data on the observables ( O exp i ). In order to determine the agreement with data, one builds a likelihood comprising the probability distributions of experimental data, evaluated at the theoretical predictions. Schematically, we multiply the probability distribution functions (pdf) provided by the experimental data in which the theoretical predictions depend on a set of given input parameters p, all associated with additional sources of uncertainty. Maximising this likelihood function then leads to the maximum likelihood estimator -i.e. "best-fit point" -as the point of highest probability. In practice, one is only interested in a subset of the theoretical input parameters, or fit parameters ( θ), leaving the remaining input parameters as nuisance parameters ( ξ) to be "integrated out". To do this, one in general follows either the Bayesian or the Frequentist approach, both computationally very expensive. Another much faster approach which is used throughout this work is a gaussian approximation of the likelihood, which can be written as In the above, O exp are the central values of the observables as measured by experiments, O th ( θ) the central values of the theoretical predictions with respect to the nuisance parameters (but dependent on the fit parameters θ i ), C exp the covariance matrix of the measurements of all included observables and C th the covariance matrix of the predictions of all included observables. The theoretical covariance matrix now contains all theoretical uncertainties of the observables (and their correlations) and is obtained by randomly sampling the nuisance parameters according to their probability distributions. Note that in this way the nuisance parameters ξ are "effectively integrated out" and the likelihood function to be optimised only depends on the parameters of interest, θ. This approach was first employed in [262]. The experimental covariance matrix is estimated by first sampling all experimental probability distributions (with a sample size of 10 6 random values), including the effects of correlations among them. In a second step, the mean values and the combined covariance matrix are estimated from the random samples. This however leads to an incorrect inclusion of strict upper limits, for instance a halfnormal distribution, since mean values of samples drawn from a half-normal distribution (or related distributions) do not correspond to the true central values, which are 0. To circumvent this problem, all observables that only have experimental upper bounds are not included in Eq. (A.2). Their likelihood is evaluated using their specific probability distributions (as provided by the experiments), at the expense of neglecting theoretical uncertainties. The probability distributions are then subsequently added to the global likelihood.
To take into account the theoretical uncertainties and correlations we use a similar Monte-Carlo method -all input parameters are randomly sampled (N MC SM = 10 4 ) according to their probability distributions. Then all observables are computed for each sample, to estimate the theoretical covariance matrix, which then also includes the theoretical correlations between observables.
The resulting approximate log-likelihood (or χ 2 ) is then minimised using the MIGRAD algorithm implemented in the minuit [263] library. For the fits of the Wilson coefficients we compute the asymmetric errors with the MINOS algorithm. For leptoquark fits this however requires excessively large computation times. Therefore, we sample the likelihoods depending on the leptoquark couplings employing MCMC-simulations using the emcee python package [261]. This results in posterior distributions of the couplings and observables of interest. The quoted 90% ranges are derived from the histograms of the posterior distributions. Here we take symmetric intervals between the 5 th and 95 th percentiles, while predicted upper limits (denoted as dashed lines) correspond to the 90 th percentile.

B Observables and data taken into account leading to the fits
In this appendix we list the observables taken into account in the different fit set-ups, as well as the datasets used for the fits. The observables (and datasets) are sorted according to the different hadronic and leptonic systems. Relevant expressions for the computation of the observables can be found in Appendix C.

B.1 Observables from b → s transitions
Leading to the fits of Sections 2.2 and 3.2, we include a large number of different binned and unbinned observables into the respective likelihoods. These play a crucial rôle in efficiently constraining the b → s transition FCNC operators and subsequently the leptoquark couplings involved.

Binned observables in b → s
We take into account all available data for the angular observables in the optimised basis [201]. Depending on the experiment providing the data, the (sub)sets of observables and bins vary. The datasets for the angular observables taken into account is summarised in Table B.1, whereas the data on the differential branching fractions is shown in Table B.2. We notice that in all cases we neglect the bin between 6 and 8 GeV 2 as, due to the cc resonances, QCD factorisation is no longer a good approximation in this region [264]. Furthermore, we do not take into account the bin [0.1, 0.98] GeV 2 : the different form factor treatments in flavio [199] and Ref. [201] lead to significant discrepancies in the associated theoretical uncertainties in this bin, while for all other bins there is a good agreement. Moreover, in the region of large hadronic recoil, we always take into account the narrow bins, whereas at low hadronic recoil we average over the kinematic region above the resonances.
In addition to the binned observables in b → sµµ, we also include the b → s LFUV observables into the likelihoods. The bins and datasets of the ratios of (differential) branching fractions R K ( * ) , as well as differences of angular observables between electrons and muons in the final state, which have either been found to be consistent with the SM, or are yet to be observed. Meson decay modes without a hadron in the final state suffer from significantly smaller hadronic uncertainties, since QCD corrections can be absorbed into a redefinition of the decay constant, and all QED and electroweak corrections remain fully perturbative. Consequently, these decays provide very clean probes for NP effects especially in C ( ) 7,10 , but also in C ( ) S,P Wilson Coefficients. A recent LHCb analysis [186] of B (s) → ee yields upper bounds at the O(10 −9 ) level. For B (s) → µµ, the situation is more complicated, since the decays are always measured in correlation to each other. While the decay B s → µµ has been observed and measured by several experiments [29,30,[202][203][204][205], as of today only upper limits on the decay B 0 → µµ are available (at the 10 −10 level), due to insufficient statistics. In order to avoid losing important correlations in the measurements, we use the 2-dimensional likelihoods (including negative values for BR(B 0 → µµ)) and sample them to obtain a naïve combination, following the prescription of Ref. [32,206].
Other observables To constrain contributions to C ( ) bsγ 7 in the dipole operator, we also include the branching fractions BR(B → K * γ) [207], BR(B → X s γ) [208] and BR(B s → φγ) [209,210]. Notice that all these observables correspond to the full branching fractions, implying that they are calculated and measured over the full kinematic region.

B.2 Charged current B-decays
Observables in b → c ν First and foremost we include the very relevant LFUV ratios R τ D ( * ) , commonly denoted R D ( * ) , into the global likelihoods. Analogously, a ratio comparing muons and electrons in the final state (R µe D * ) can be defined, which shows excellent agreement with the SM [211,212]. For R τ D * we use the uncorrelated measurements by LHCb [9,213] and Belle [10], whereas for R τ D there are several measurements, obtained by BaBar [4] and Belle [6,10,11], always in correlation with R τ D * . Numerous other observables are taken into account in addition to the anomalous ratios R D ( * ) . The extensive array of experimental data (in binned branching fractions of the decay B → D ( * ) ν) used in our fits is presented in Table B.4. Furthermore, we include the unbinned branching fractions BR(B + → D ( * ) µν), BR(B + → D ( * ) eν) [214,215] and the inclusive branching fraction BR(B → X c eν) [216,217].
Other charged current B-decays In addition to charged current b → c ν decays, we also include certain b → u ν decays to obtain further constraints on the leptoquark couplings to the first quark generation. These can be found in Table B.5.

B.3 Strange, charm and τ -lepton decays
The above listed data mostly allows to constrain combinations of second and third generation quark leptoquark couplings (to all leptons). To achieve more precise constraints for the second and first generation quarks, we further include numerous decays of strange and charm flavoured mesons. Since the light mesons cannot decay into τ -leptons, we also use data on SM allowed τ -lepton decays, as a complementary source of information.  [199].
Binned charm decays In addition to the precise measurements of the full branching fractions of several charmed meson decay modes, there are also precise measurements of the q 2 distributions for several charged current decay modes in semileptonic charm decays with an electron in the final state. The datasets used are presented in Table B Unbinned observables Besides the binned semileptonic charm decays, we also include the full branching fractions for charged current leptonic and semileptonic charm decays, charged and neutral current decays of strange flavoured mesons, and charged current semileptonic τ -lepton decays. The charged current decays are listed in Table B.7 and the neutral current ones in Table B.8. Table B.7: Data on charged current charm and strange flavoured meson decays. The SM predictions are obtained using flavio [199].

B.4 Belle II Observables
As discussed in Section 4.2, we use specific fit set-ups which allow for an extrapolation of the current situation into the near future. The future sensitivities, taken into account as data, are listed in Table B.9; these always correspond to the full anticipated luminosity of 50 ab −1 .

C Vector leptoquark contributions to leptonic and mesonic flavour observables
New physics models aiming at addressing the LFUV hints in B-meson decays typically give rise to new contributions to several flavour observables depending on the new flavour structure; these include contributions to various flavour conserving and flavour violating leptonic and semileptonic mesonic decay modes, as well as cLFV processes. In particular, the vector leptoquark scenario can already  [15,22] (1.20 ± 0.12) × 10 −7 [142] < 2 × 10 −5 Table B.9: Observables for which Belle II will improve on current experimental sensitivities. The SM predictions are obtained using flavio [199], unless otherwise stated.
contribute to some of these observables at tree level, while others receive leading contributions at the one-loop order. In this appendix we collect information allowing to estimate the vector leptoquark contribution to several of the above mentioned processes.

C.1 Leptonic and semileptonic meson decays
Here we summarise the different vector leptoquark contributions to leptonic and semileptonic meson decays which arise at tree-level, and to modes with final state neutrinos (whose new contributions arise at one-loop level). We do not include neutral meson oscillations which arise at one-loop level and typically provide much weaker constraints if, apart from the leptoquarks, only SM fields are considered. However, we notice that this may no longer hold in the presence of additional heavy fermionic states (which might be present in a UV-complete model, as for example heavy vector-like leptons); in that case, the contributions could be sizeable so that neutral meson oscillations can then lead to important constraints, as discussed in [96,246].
Vector leptoquarks can induce new contributions to purely leptonic decays of pseudoscalar mesons, leading to important constraints on the flavour structure of V 1 couplings. Here, we provide a brief summary of the formalism for the computation of the P → − + rates. Following the standard decomposition of the hadronic matrix element [283] 0 |d j γ µ γ 5 d i | P (p) = i p µ f P , where f P corresponds to the P meson decay constant, the branching fraction can be expressed as where the λ(a, b, c) is the standard Källén-function, defined as λ(a, b, c) Note that for a lepton flavour conserving decay mode, e.g. B s → µµ, one must include the SM contribution and the relevant RG running effects. Since the vector leptoquarks contribute to the leptonic pseudoscalar meson decays at the tree level, such processes can provide important and very stringent constraints on the vector leptoquark couplings.
The semileptonic decays of pseudoscalar mesons can also be the source of significant constraints on the vector leptoquark couplings. To evaluate the differential branching fractions for these modes, we parametrise the hadronic matrix elements following the standard convention as where the momentum transfer lies in the range (m + m ) 2 ≤ q 2 ≤ (M P − M P ) 2 . For the evaluation of the form factors we closely follow the prescription of [161]. The final differential branching fraction for the decay P → P − + can be expressed in the form d BR(P → P − + ) dq 2 = |N P (q 2 )| 2 × ϕ 7 (q 2 ) |C 7 + C 7 | 2 + ϕ 9 (q 2 ) |C 9 + C 9 | 2 + ϕ 10 (q 2 ) |C 10 + C 10 | 2 + ϕ S (q 2 ) |C S + C S | 2 + ϕ P (q 2 ) |C P + C P | 2 + ϕ 79 (q 2 ) Re (C 7 + C 7 ) (C 9 + C 9 ) * + ϕ 9S (q 2 ) Re (C 9 + C 9 ) (C S + C S ) * + ϕ 10P (q 2 ) Re (C 10 + C 10 ) (C P + C P ) * , and the normalisation factor is given by The vector leptoquark can also contribute to s → dνν and b → sνν transitions at one-loop level. The |∆S| = 1 rare decays K + (K L ) → π + (π 0 ) ν ν and B → K ( * ) ν ν correspond to the quark level transition d j → d i ν ν , which can be described by the short-distance effective Hamiltonian [284][285][286] where i, j corresponds to the down-type quark content of the final and initial state mesons, respectively. For vector leptoquarks, the one loop contributions are a priori divergent; consequently, the corresponding would-be Goldstone modes must be consistently included to obtain the correct result.
Following the prescription of [185], the coefficient C ij L,f a for d a → d fνi ν j , due to V 1 leptoquark exchange is given by where M W and m t respectively correspond to the masses of the W boson and top quark. The neutral and charged kaon decay branching fractions can then be obtained by [276,287] BR(K ± → π ± νν) = 1 3 (1 + ∆ EM ) η ± × (C.14) Here, λ corresponds to the standard Wolfenstein parametrisation (i.e. the Cabibbo angle), λ c = V * cs V cd and λ t = V * ts V td . The decay width for B → K ( * ) νν has been derived in [284], leading to C SM,f i L,sb ≈ −1.47/s 2 W δ f i , which can be used to normalise the branching ratios as

C.2 Charged lepton flavour violating decays
Charged lepton flavour violating observables, such as radiative decays i → j γ, three-body decays i → 3 j , as well as neutrinoless µ − e conversion in nuclei, can lead to important constraints on the vector leptoquark couplings, due to the non-universal couplings to different flavours of SM charged leptons. We recall that while i → j γ and i → 3 j decays can be induced at one-loop level by the vector leptoquark, µ − e conversion in nuclei can occur at tree-level. Here also, the one-loop dipole and anapole contributions from the exchange of a vector leptoquark are a priori divergent and to obtain a finite result the would-be Goldstone boson degree of freedom (degenerate in mass with vector leptoquark) must be included. After symmetry breaking, the latter degree of freedom is subsequently absorbed by the massive vector leptoquark.
(C.20) 9 As discussed in Section 3.2, we recall that in the current study we work under the assumption that K ij R 0.

C.2.2 Three body decays →
Vector leptoquarks can induce three-body cLFV decays → at the loop level, through photon penguins (dipole and off-shell "anapole"), Z penguins and box diagrams. The effective Lagrangian relevant for these decays can be expressed as [289,290] L → eff = L i → j γ eff − 4 G F √ 2 g 1 (¯ P L )(¯ P L ) + g 2 (¯ P R )(¯ P R ) + + g 3 (¯ γ µ P R )(¯ γ µ P R ) + g 4 (¯ γ µ P L )(¯ γ µ P L ) + + g 5 (¯ γ µ P R )(¯ γ µ P L ) + g 6 (¯ γ µ P L )(¯ γ µ P R ) + H.c. , (C. 21) where the photonic dipole part, cf. Eq. (C.16), with the corresponding Wilson coefficients C i j L(R) have already been discussed in the previous subsection; the off-shell anapole photon penguins, Z penguins and box diagrams contribute to g 3 , g 4 , g 5 and g 6 coefficients. For our numerical analysis we only include the log-enhanced photonic anapole contributions 10 in addition to the dipole ones. In the absence of right-handed couplings of the vector leptoquark, the only non-vanishing coefficients are g 4 = g 6 given by where Q f = Q denotes the charge (in units of e) of the fermion pair attached to the end of the off-shell photon and with the loop function f a (x) given by In the above, N c denotes the number of colours of the internal fermion and x i = m 2 d i /m 2 V 1 . As an example, in the case of µ → 3e decays, the branching ratio can be written as [289,290] BR(µ → eee) = 2 |g 3 | 2 + |g 4 | 2 + |g 5 | 2 + |g 6 | 2 + +8 e Re C µe R (2g * 4 + g * 6 ) + C µe L (2g * 3 + g similar expressions for the other cLFV 3-body decay modes can be obtained in a straightforward manner.

C.2.3 Neutrinoless µ − e conversion
Neutrinoless µ − e conversion can be induced by the vector leptoquark V 1 at tree level, in addition to the one-loop contributions through dipole and anapole photon penguins. Therefore, µ − e conversion provides very stringent limits on the vector leptoquark couplings to the first two generations of SM charged leptons. The general contribution to the neutrinoless µ−e conversion due to vector leptoquark can be written as [197] Γ(µ − e, N) = 2 G 2 F C µe * R m µ D + 2 g (u) LV V (p) + g (u) LV + 2 g where the photonic dipole Wilson coefficients C i j L(R) can be found in Eq. (C. 18) and (C.19); the other non-vanishing Wilson coefficients, relevant for vector leptoquark exchange, are given by Here, Q d = − 1 3 , Q u = 2 3 , and the values for the overlap integrals (D, V, S) can be found for instance in [292]. The relevant scalar coefficients G (d i ,N ) S are given in [293].