Flavorful leptoquarks at the LHC and beyond: Spin 1

Evidence for electron-muon universality violation that has been revealed in $b\to s \ell\ell$ transitions in the observables $R_{K,K^*}$ by the LHCb Collaboration can be explained with spin-1 leptoquarks in $SU(2)_L$ singlet $V_1$ or triplet $V_3$ representations in the ${\cal{O}}(1-10)$ TeV range. We explore the sensitivity of the high luminosity LHC (HL-LHC) and future proton-proton colliders to $V_1$ and $V_3$ in the parameter space connected to $R_{K,K^*}$-data. We consider pair production and single production in association with muons in different flavor benchmarks. Reinterpreting a recent ATLAS search for scalar leptoquarks decaying to $b \mu$ and $j \mu$, we extract improved limits for the leptoquark masses: For gauge boson-type leptoquarks ($\kappa=1$) we obtain $M_{V_1}>1.9$ TeV, $M_{V_1}>1.9$ TeV, and $M_{V_1}>1.7$ TeV for leptoquarks decaying predominantly according to hierarchical, flipped and democratic quark flavor structure, respectively. Future sensitivity projections based on extrapolations of existing ATLAS and CMS searches are worked out. We find that for $\kappa=1$ the mass reach for pair (single) production of $V_1$ can be up to 3 TeV (2.1 TeV) at the HL-LHC and up to 15 TeV (19.9 TeV) at the FCC-hh with $\sqrt{s}=100$ TeV and $20 \, \mbox{ab}^{-1}$. The mass limits and reach for the triplet $V_3$ are similar or higher, depending on flavor. While there is the exciting possibility that leptoquarks addressing the $R_{K,K^*}$-anomalies are observed at the LHC, to fully cover the parameter space $pp$-collisions beyond the LHC-energies are needed.


I. INTRODUCTION
The B-decay observables R K and R K * [1] probe electron-muon universality violating new physics (NP) in flavor changing neutral current b → s transitions. Recent measurements by the LHCb Collaboration [2,3] returned values that are both about (15 − 25)% lower than the lepton-universality limit R K = R K * = 1, each with statistical significance up to 2.5 σ. Very recently, the LHCb Collaboration presented a new update of R K [4] with the central value equal to their previous result [2] but with reduced uncertainty, namely: R LHCb(2021) K = 0.846 +0.044 −0.041 , revealing a 3.1 σ deviation [4] from the SM expectation. Naive combination of R K and R K * measurements shows a deviation from the standard model above 4σ. Even more precise measurements are therefore mandatory to clarify these anomalies. In addition, correlations with other observables are to be expected if lepton universality is actually violated in semileptonic B-decays. One important cross check is to look for lepton universality breaking in other rare decay channels [5], c.f. recent tests in baryonic decay mode Λ b → pK [6], which intriguingly point in the same direction as the ones in B → K ( * ) -a suppression of dimuon modes relative to dielectron ones -although with lesser significance. The anomalies also suggest to undertake direct searches for the corresponding NP particles at the Large Hadron Collider (LHC) and future high-energy proton-proton (pp) colliders [7][8][9][10].
An interesting class of models that can naturally accommodate the measured values of R K,K * involves a heavy leptoquark that couples to b -and s currents, thus contributing to b → s transitions at tree-level, with masses at scales of up to a few tens of TeV. Three suitable representations are singled out using both R K and R K * , that is, the SU (2) L -triplet scalar S 3 (3, 3, 1/3), and two spin-1 multiplets, the SU (2) L -singlet V 1 (3, 1, 2/3) and the triplet V 3 (3, 3, 2/3) [5,[11][12][13][14][15]. Specifically, R K and R K * constrain the leptoquark's coupling to strange quarks and leptons , λ s , and to b-quarks and leptons, λ b , divided by the square of the leptoquark mass M V as where we included R K [4] combined with R K * [3] in the q 2 -bin (1.1, 6) GeV 2 , see Refs. [5,14] for details 1 . One verifies that i) this indeed requires the leptoquark to couple differently to muons than to electrons, that is, lepton nonuniversality, ii) that the data point to a collider mass scale, and iii) that further input is required to extract values of the individual leptoquark couplings, and therefore, the mass. The latter can be in reach of the LHC if the couplings are sufficiently small. In general, given sufficient energy, a collider study can not only discover leptoquarks consistent with (1), but also determine the leptoquark couplings and mass.
In this paper we focus on the collider signatures of spin-1 (vector) leptoquarks, as a sequel to [16] on the scalar leptoquark S 3 . As in our previous works we assume that the dominant lepton species involved in (1) are muons. This choice is pragmatic, as both lepton species could couple to NP and sizeable couplings of leptoquarks to electrons are presently not excluded. However, choosing muons over electrons as the main contributors to (1) is consistent with the global b → s fits 2 , which also suggest NP in b → sµµ angular distributions. We stress that dielectron decays b → see also deserve precise experimental treatment in the future [17,18].
To quantitatively study the sensitivity to leptoquarks at the LHC and beyond, the leptoquark couplings to bottomand to strange quarks have to be given individually, not only their product as in (1). This is obvious for single production, which feeds on the corresponding parton distribution functions (pdfs) in the proton, see Fig. 1, but matters also for pair production [19], since the flavor patterns dictate the signature from leptoquark decay. Flavor symmetries [20] can provide such an input to leptoquark couplings [21], implying hierarchical pattern with λ s /λ b proportional to the strange over the b-quark mass. To explore more general settings we employ in addition a flipped benchmark pattern, with inverted hierarchy, and a democratic one. Note that flavor non-diagonal couplings to leptons and quarks are required to explain the anomalies (1), as is taken into account in recent leptoquark searches in pair production at ATLAS [22,23]; see [24][25][26][27] for other recent single and pair production searches and [28] for a review on leptoquark phenomenology. This paper is organized as follows: In Sec II we introduce the relevant interaction terms for V 1 and V 3 , and introduce three flavor scenarios for the relative size of the couplings to second and third quark generations. In Sec. III we discuss single-, pair-and resonant production mechanisms in the final-state channels involving muons, strange-or bottom quarks based on (1) and the flavor benchmarks. We work out present mass limits and determine the sensitivity at the HL-LHC and future pp-colliders. Going beyond the study [16] of the S 3 leptoquark, here the projections are based on extrapolations of existing LHC searches. We conclude in Sec. IV. In the appendix we give an approximate analytical argument for comparing V 1 and V 3 single production cross sections.

II. MODEL SETUP
We briefly review the vector leptoquarks V 1 and V 3 and their interactions with standard model particles in Section II A. In Section II B we discuss predictive flavor patterns in the context of the R K ( * ) anomalies.
We start by recalling the new interaction terms that need to be added to the SM Lagrangian. The couplings of V 1 to leptons and quarks are 1 Current data results in the ratio [5] X K * ≡ R K * /R K = 0.84 ± 0.13, consistent with X K * = 1 at ∼ 1.3 σ, which indicates the possibility of a small admixture from right-handed currents. These could stem from scalarS 2 (3, 2, 1/6) or vector leptoquarks V 2 (3, 2, −5/6). In a two-dimensional fit with both left-and right-handed currents the right-hand side of Eq. (1) reads (−1 ± 0.23)/(37 TeV) 2 . 2 It is possible that both muon-and electron channels are affected by contributions of opposite sign [53]. respectively. Here, T a denote the generators of QCD in the fundamental representation, normalized to tr(T a T b ) = δ ab /2, and V 3 , W µν are the three-component vectors in spin-1 representation of SU (2) L . In addition to the terms with covariant derivatives D µ renormalizable, gauge invariant interactions with the gluon field strength tensor, parametrized by the coupling κ, exist (see e.g. [29,30]). Other renormalizable couplings of the leptoquark bilinear to weak boson field strengths W µν and B µν are irrelevant for the present study. The value of κ depends on the ultraviolet completion of the model, e.g. κ = 1 in a Yang-Mills case in which the vector leptoquark is the gauge boson of a non-abelian gauge group. We choose the benchmark values κ = 0 and κ = 1 throughout the paper, see Sec. III A for a brief discussion of the impact of κ on leptoquark production.
We assume that only couplings to quark and lepton doublets are present and hence neglect couplings to singlet fermions in Eq. (2). This feature is not generic across possible UV completions and requires some model building, e.g. it does not hold in the minimal Pati-Salam model. Some of the proposed models in which V 1 is a gauge boson also include new vector-like fermions which render the coupling matrix λQ L non-unitary. Several models that can accommodate such a choice for V 1 have been proposed in the literature. For a selection of references studying the V 1 leptoquark in the context of the b → s transitions we refer the reader to Refs. [13,.

B. Three flavor benchmarks
The measured values of R K [4] and R K * [3] can be accommodated with the combination of couplings and the leptoquark mass given in (1). It is apparent that an additional constraint on the leptoquark's parameter space with regards to collider searches is required. For instance, λ bµ λ sµ ∼ 1 points to a mass scale around 40 TeV, outside of the search range of any presently planned collider, whereas weaker couplings λ bµ λ sµ ∼ 10 −3 −10 −2 imply lower leptoquark mass, in reach of the LHC. Flavor symmetries, which explain the observed pattern of standard model masses and mixings, do provide naturally requisite suppression mechanisms [13]. These symmetries determine the ratio between the leptoquark couplings to b-and s quarks. We employ these constructions when defining flavor benchmark pattern. In addition, the B s -B s mass difference, to which V 1,3 contribute at 1-loop, combined with R K,K * , impose upper bounds of around 45 TeV and 20 TeV on the masses of the V 1 and V 3 leptoquarks, respectively [14]. More recent analysis of B s -B s mixing constraints on NP in R K,K * suggest even lower upper mass limits Ref. [54]. The bound from the loop-induced B s -B s mixing turns out to be dependent on the specific completion of the vector leptoquark model into a renormalizable theory at high energies, contrary to models with scalar leptoquarks. A comprehensive fit to available leptonic and semileptonic b → sµµ data sharply supports the C 9 = −C 10 solution [55] that corresponds to our models [14] and shows consistency with R K,K * and other muon specific observables. Our scenarios also induce lepton universality violation in the charged current semileptonic B decays within the ratio R 1/(5.4 TeV) 2 which is safely satisfied by our relation (1). Furthermore, the constraints from perturbative unitarity place an upper bound of about 80 TeV on the scale of effective operator relevant for our present setup [57].
In the following we consider three benchmark scenarios with coupling textures that couple the vector leptoquark predominantly to the second lepton generation.
Hierarchical scenario The first texture we consider is the same as in Ref. [16] based on flavor models discussed in Ref. [13], where we assume that the hierarchies found in the standard model masses and mixings are also present in the leptoquark couplings. This is the case in simple flavor models based on the Froggatt-Nielsen-Mechanism [20] which induces the hierarchies between the different quark generations, where ∼ 0.2 is of the order of the Wolfenstein parameter, i.e. the sine of the Cabibbo angle. Specifically, we employ where we assume contributions to the first quark generation to be suppressed strongly enough to not violate any existing bounds from data on µ-e conversion or rare kaon decays. Entries marked with " * " arise only through higher order corrections within the models from Ref. [13]. The parametric suppression of the individual quark generations is preserved by Cabibbo-Kobayashi-Maskawa (CKM) rotations. As neutrinos are reconstructed inclusively at collider experiments flavor rotations in the lepton sector do not affect such observables.
Allowing for an additional O(1) factor in the ratio between λ sµ and λ bµ couplings, taken within 1/3 and 3, the central value on the right handed side of Eq. (1) implies These additional flavor model uncertainties dominate over the experimental ones in (1).
Flipped scenario As a second scenario, we consider the inverted form of the previous texture, that is: This yields the same effect in b → s + − transitions as the hierarchical pattern while enhancing the single production cross section due to larger pdf of the strange quark. We obtain the same coupling range for λ 0 as in the hierarchical scenario given in Eq. (11).
Note that this pattern has a weaker foundation in flavor models, and if it is introduced in the interaction basis the CKM rotations can induce contributions to first generation quarks at order .
Democratic scenario Lastly, we consider a texture where the couplings to the second and third quark generation are of equal size: Taking into account the aforementioned O(1) factor and Eq. (1) imply Each scenario contains four parameters, the mass, the parameter κ and the dominant couplings λ bµ and λ sµ . The measurements of the single-or pair-production cross section, the corresponding branching fractions and the resonance width, together with the reconstruction of the mass peak, would suffice to determine all four parameters. Note that b-tagging would be necessary for such an analysis.
We stress that flavor models which isolate a single species of leptons are straightforward to obtain using techniques from neutrino model building [21], however, more general settings can arise and are viable, too. Allowing for significant entries " * " in (10), (12), (13) would open up further leptoquark decay modes and search channels, which would reduce branching ratios in the signal channels studied here. Negligible entries " * " correspond therefore to the most favorable situation for an observation in the muon channel.
Leptoquarks V 1,3 and S 3 contribute also at tree level to charged current b → c ν-decays [12,13,[32][33][34][35][36][37][38][39][40][41][42][43][44]58]. In particular V 1 , which evades dominant constraints from b → sνν on third generation lepton couplings, has received interest as a possible resolution of the anomalies in the b → cτ ν observables R D ( * ) . Note that (1) points to a NP contribution to a loop-induced process in the standard model, hence a corresponding NP effect in tree level charged currents would naturally be O(1/(4π) 2 · λ qτ /λ qµ )-suppressed. An effect of the same order of magnitude in the charged current as in the neutral current, which is about ten percent, would therefore require substantial hierarchy λ qτ ∼ 10 2 λ qµ , which is unsupported by flavor models and points to strong couplings to τ 's. As a study in concrete, full flavor models [13] is beyond the scope of this work, we also do not consider links with the b → cτ ν observables.

III. COLLIDER PHENOMENOLOGY
In this section we study vector leptoquark production in pp-collisions and decays of leptoquarks. Basics are given in Sec. III A. We work out bounds on the masses of vector leptoquarks using available search results for pair-production of leptoquarks from ATLAS [22] (Sec. III B), and cross sections for future pp colliders in Sec. III C. We consider three setups corresponding to center-of-mass energies √ s: 14 TeV (LHC run 3), 27 TeV (HE-LHC), and 100 TeV (FCChh) [8] with target integrated luminosities of L = 3 ab −1 , 15 ab −1 and 20 ab −1 , respectively. In Sec. III D we also briefly discuss resonant production. We analyze the mass reach of future pp colliders by extrapolating current limits on cross sections to higher center-of-mass energies and luminosities in Sec. III E.

A. Leptoquark production and decay
We consider three dominant mechanisms of leptoquark production at pp colliders: pair production, single production in association with a lepton and resonant-production induced by quark-lepton fusion, shown in Fig. 1. The structure of the final state signatures of all three classes of processes is determined by the flavor structure of the leptoquark couplings.
The flavor scenarios (10), (12), (13) can be distinguished experimentally by different patterns of the final states in two-body decays of the leptoquarks. In the hierarchical scenario (10), the dominant leptoquark decay modes are for the singlet and for the triplet. The bµ + and tν final states of the V +2/3 1,3 leptoquarks are related by SU (2) L symmetry such that their branching fractions are approximately equal.
In the flipped scenario (12) the leading signatures involve charm and strange quarks for the singlet and for the triplet. In the democratic scenario (13) all of the above modes arise and final states with both light and heavy quarks are relevant. We recall that we allow in (13) a mild hierarchy between the two couplings which can have a strong impact on the relative size of the different final state branching ratios as B ∼ |λQ | 2 . Approximate branching ratios for the benchmark patterns (10), (12), (13) are given in Tables I and II. bµ + tν sµ + cν  In Fig. 2 we show the pair-and single production cross sections as functions of κ for the example of the HL-LHC √ s = 14 TeV and M V1 = 3 TeV. The cross sections exhibit minima for κ in the vicinity of 0 or −1, see also [29]. The shapes vary mildly with the variation of the leptoquark mass in a range suitable for a 14 TeV collider.

B. Current mass bounds
Leptoquark-based explanations of the deviations found in B-physics motivated the recent search for pair-produced scalar leptoquarks with 139 fb −1 of data from 13 TeV p-p collisions [23]. Mass limits for scalar leptoquarks decaying (7) of the single-(red, solid) and pair production cross section (green, dashed) for V1. We fix √ s = 14 TeV and MV 1 = 3 TeV. For the single production cross section we employ the hierarchical scenario. Analogous results are obtained for other choices of the parameters and flavor benchmarks.
dominantly to top and electron (muon), obtained from this search, are 1470 GeV (1480 GeV). Such final states appear in decays of scalar (vector) SU (2) L triplet leptoquark S 3 (V 3 ). Previously derived collider bounds on vector leptoquarks from pair production searches are M V1 > 1.3 TeV for the dominant decays to τ b, and M V1 > 1.7 TeV in the µb channel [59], for κ = 1.
We evaluate the current mass limits for the V 1 and V 3 leptoquarks using the limits on the cross sections found in the ATLAS Collaboration search [22] for pair production of scalar leptoquarks. This search was performed using the data collected in the 13 TeV LHC run with a luminosity of L = 139 fb −1 . For the hierarchical scenario we use the limits obtained in the (bµ , bµ)-channel while for the flipped scenario we use the (qµ , qµ)-channel. The role of q in the latter channel is played in our V 1 model by the strange quark, and by charm for V 3 . The bound in the democratic scenario is obtained from the (bµ , bµ)-channel. The limits on the cross section in the (bµ, bµ)-, (qµ, qµ)-and (cµ, cµ)-channels (bottom plot) are shown in Fig. 3 (plots to the left) -the comparison to the theoretical cross sections for leptoquark pair production and subsequent decay in corresponding final-state channels determines the mass limits.
We now spell out the obtained limits for the V 1 -leptoquark model. For both the hierarchical and the flipped scenarios the limits are 1.7(1) TeV and 2.0(1) TeV for κ = 0 and κ = 1, respectively. The bounds are the same for both of these scenarios because the current experimental limits for the cross sections to (bµ, bµ)-and (sµ, sµ) final states nearly coincide in the region of large leptoquark masses. The limits are somewhat weaker for the case of democratic scenario and read 1.5(1) TeV and 1.8(1) TeV for κ = 0 and κ = 1, due to smaller individual branching fractions into bµ-and sµ-pairs, see Tab I. For the democratic scenario, in which the leptoquark couplings to b− and s quarks are approximately equal and both of the above mentioned limits apply, we used (bµ, bµ)-channel, since the corresponding experimental limit is currently somewhat more strict in the most of the explored mass range than the one for (qµ, qµ).
In the case of the model with the V 3 leptoquark in the hierarchical scenario the limit is the same as for the case of The mass limits depend on the value of κ, as shown in Fig. 3 (plots to the right). Note that the cross sections have a minimum within κ ∈ (−2, 2), which corresponds to the weakest bound on the mass. These are the same for V 1 and V 3 in the hierarchical scenario, M V1,3 > 1.6 TeV for κ = −0.3, and equal to one for V 1 with flipped and V 3 with democratic flavor structure. Note, however, that the absolute minimum for the bound on the V 1 mass is obtained in the democratic scenario and reads M V1 > 1.4 TeV for κ = −0.3. The corresponding weakest bound in the case of V 3 is found in the case of hierarchical scenario, given above.
Constraints on the parameter space of leptoquark models can also be obtained from Drell-Yan processes, to which the leptoquarks contribute via t-channel exchange [60,61]. For our present setup, the corresponding di-muon channel is relevant. The authors of Refs. [44,61] performed the recast of the ATLAS collaboration measurement [62] in the dimuon channel at 13 TeV with 36 fb −1 of data, while the authors of Ref. [50] used the result by the CMS collaboration [63] corresponding to the same center-of-mass energy and luminosity. The results of these studies imply  [22]. Top row: V1 in the hierarchical and flipped flavor scenarios, which equals V3 in the hierarchical and democratic flavor scenarios, middle row: V1 in the democratic scenario, and bottom row: V3 in the flipped scenario. Left: Sensitivity to V1, V3-pair production, assuming dominant decays to µb, µq and µc in purple, blue and black, respectively; q denotes quarks lighter than the charm quark, The green bands indicate the theory prediction including the pdf-and scale uncertainties, for κ = 0 (solid) and κ = 1 (dashed). Right: Mass bounds for the leptoquarks V1, V3 as the function of parameter κ (7). The boundary of the excluded region is represented by the band whose width results from the pdf-and scale uncertainties. that the parameter space relevant for our model does not receive constraints at present. This situation could change as more data is collected in the future [61].
We find that the cross sections for pair production of the V 1 leptoquark typically turn out several times larger than those for any specific weak-isospin component of the scalar leptoquark S 3 studied in Ref. [16]. This is not surprising, given that the vector leptoquark involves three helicity states. Thus, for given values of the branching fractions in the specific lepton-quark channels, taken to be equal for the scalar and the vector leptoquark, the corresponding search limits for the scalar turn out weaker, in accord with what was previously noted in Ref. [44]. However, we note that the pattern of the branching fractions into the final state lepton-quark pairs is guided by the details of the flavor structure of the Yukawa couplings to the fermions, and the SU (2) L structure of a leptoquark representation. For example, for M LQ = 3 TeV and √ s = 14 TeV, we have σ(pp → S → bµ) 1/2, which lowers the cross section of the vector pair-production in the (bµ, bµ) channel by factor 1/4.

C. Single and pair production cross sections
We evaluate the leading order cross sections for the single production of V 1 in association with a muon, represented by the resonant diagrams (a) and (b) in Fig. 1, as functions of the leptoquark mass. The results are displayed by the red bands using κ = 0 in Fig. 4 for the case of V 1 , and in Fig. 5 for V 3 . The results for κ = 1 are displayed by hatched bands. The leading-order cross sections for pair production and subsequent resonant decays are represented by the solid (hatched) light green bands for κ = 0 (κ = 1).
For each flavor scenario we assume that the parameters of the leptoquark model satisfy Eq. (1). The band widths originate from Eqs. (11), (14). We note that there are also non-resonant diagrams contributing to (q , q ) final states that were not taken into account in our numerical analysis, assuming that the contributions of the resonant diagrams shown Fig. 1 (c)-(g) are well separated by appropriate kinematic cuts.
For the evaluation of the cross sections and the corresponding uncertainty bands we used Madgraph [64] with the UFO [65] output of the leptoquark models that we implemented using Feynrules [66]. The pdf-, and scale uncertainties are evaluated using LHAPDF [67], symmetrized and combined in quadrature. Our Feynrules implementations of the V 1 and V 3 models are attached to this paper as ancillary files. We checked the consistency of our Feynrules implementations with the corresponding implementations from Ref. [68]. We used the software package Feyncalc [69,70] for several cross-checks.
Pair production is predominantly induced by the QCD-initiated processes and is essentially independent of the flavor structure. The latter determines the branching fractions into various final state channels, see Tab. I. As an exception to this, the large contribution of diagram (e) shown in Fig. 1 becomes noticeable for the large-mass region within the flipped scenario, see the last plot of the second row in Fig. 4. In this case, the R K,K * -condition (1) for flipped flavor hierarchy, forces the large values for the coupling λ sµ , reaching the borders of the perturbativity. This as well as the larger uncertainty band in this plot can be understood from (11) for large values of M V .
The magnitude of the single production cross section induced by qg → V 2/3 1 at parton level is directly proportional to the square of the magnitude of the corresponding flavor coupling λQ . Assuming the narrow width approximation, we multiply the corresponding production cross sections by the corresponding branching fractions given in Tab. I. Since there are no available single production searches involving b quarks in the final state, we added the contributions involving jets and b-quarks which amounts to the branching fraction 1/2 for each of the three flavor scenarios. Single production for κ = 0 in association with R K ( * ) Pair production for κ = 0 Sensitivity (pair production) Single production for κ = 1 in association with R K ( * ) Pair production for κ = 1 Sensitivity (single production) , derived from the RK,K * -band in Eq. (11) for hierarchical and flipped scenarios, and Eq. (14) for the democratic scenario. Light green: pair production with final states (bµ, bµ) for the hierarchical and democratic scenarios and (qµ, qµ) for the flipped scenario. The error bands for pair production are evaluated by combining the pdf-, and scale uncertainties. Results for κ = 1 are shown in a dashed/hatched form together with the solid curves for κ = 0. The solid dark red and the dark green curves depict the projected experimental sensitivity for single and pair production, respectively. As the starting curves for these extrapolations, we used the results of the measurements by the CMS [25] and ATLAS [22] collaborations, for the single-and pair production, respectively, see Sec. III E for details. The dot-dashed segments of the extrapolated curves for the low masses required additional smooth variation of the luminosities between the initial and the target values, following the prescription in Ref. [79]. The dotted segments for the large masses represent the smooth continuation above the final extrapolated points towards the higher masses with the constant values of the cross-section limit. Neither of these segments play a role in determining the mass reaches shown in Tab. III. Single production for κ = 0 in association with R K ( * ) Pair production for κ = 0 Sensitivity (pair production) Single production for κ = 1 in association with R K ( * ) Pair production for κ = 1 Sensitivity (single production)  Light green: pair production with final states (bµ, bµ) for the hierarchical and democratic scenarios and (cµ, cµ) for the flipped scenario. The error bands for pair production are evaluated by combining the pdf-, and scale uncertainties, see Fig. 4 and Sec. III E for the details.

D. Resonant production
Determinations of the photon distribution function inside the proton introduced in Refs. [71,72] were recently followed by the determination of the lepton pdfs in Ref. [73]. These results opened up the possibility to consider resonant leptoquark production from lepton-quark fusion in pp collisions [74], see diagram (h) in Fig. 1. Next-toleading-order QCD and QED corrections to the resonant production of scalar leptoquarks have recently become available [75]. To illustrate the expected range within the flavor scenarios (10), (12), (13), we give in Fig. 6 the resonant cross section for V 1 at a √ s = 14 TeV pp collider. The results are obtained by convolution of the leading order partonic cross section for µ (b + s) → V 1 with the LUXlep-NNPDF31 nlo as 0118 luxqed [73] pdf set that includes the leptonic pdfs. Note that the charge conjugated process is also included in the results. We parse the pdf set in Mathematica using the package ManeParse [76]. The resulting cross sections are larger than those of pair-and single production due to lesser phase-space suppression. It would be desirable to look for collider signatures of this process. Resonant vector leptoquark collider signatures motivated by the R D,D * anomalies, and the corresponding backgrounds, were recently discussed in Ref. [77]. In Fig. 7 we compare the cross sections of the resonant-and single production for V 1 -masses up to ∼ 10 TeV. The regions of the (M V1 , √ s)-plane to the right of the thick blue lines result in resonant cross section being larger than the one for single production. The current level of the lepton pdf-uncertainties does not allow for extrapolations to higher energy scales. We note in passing the absence of the triplet leptoquark component V −1/3 3 in the resonant production -its coupling to the fermion sector exclusively involves neutrinos, see Eq.(5).

E. Sensitivity projections for future colliders
In order to estimate the mass reach of the future colliders for the flavor benchmark scenarios, we extrapolate existing bounds from single-and pair production using the limit extrapolation method following Refs. [78,79]. The method assumes that the exclusion limits are determined by the numbers of background events and involves the appropriate re-scaling of the background processes with the corresponding parton luminosity functions, see [78,79] for more details. We expect the method to be less suitable for the case of leptoquarks than for e.g., the case of s-channel resonances, for which it was initially used [78], however, it should provide the correct estimate of the order of magnitude for the collider limits on the corresponding cross sections.
As the starting point for our approximation for the future sensitivity projections for single production, we employ the limits obtained by the CMS Collaboration in the √ s = 8 TeV run with L = 19.6 fb −1 [25]. The latter paper presents the limits on the resonant cross sections of the single production of the leptoquarks in association with muons in the µµj final states. Our extrapolations assume that the final b quark is not tagged and is counted as a light jet, however, we stress once again that b-tagging is required for distinguishing between the flavor scenarios, and could lead to improved limits in the case of the hierarchical and democratic flavor scenarios.
For the extrapolations of the limits on cross sections for pair production we use the search performed by the ATLAS Collaboration in Ref. [22] at 13 TeV with L = 139 fb −1 , see also Refs. [24,26] for earlier searches. We use the leading order set of pdfs provided by MSTW Collaboration [80] for the evaluation of the extrapolations. As a cross-check, we use the pdf-set NNPDF23 lo as 0130 qed [81], parsed using the package ManeParse [76]. For the hierarchical and  Mass reach in TeV for vector leptoquark single production in the hierarchical, flipped and democratic scenarios from Sec. II B and pair production, at different future colliders for κ = 0 and κ = 1. For single production we provide the mass reaches corresponding to the upper limit of the cross section band resulting from Eqs. (11), (14). In the flipped and democratic scenarios as well as for pair production we show the increased mass reaches for V3 in parentheses, for the hierarchical scenario the V1,V3 reaches are the same, see Appendix A for details.
democratic scenarios we extrapolate the limits in the (bµ , bµ)-channel, while the limits for the (qµ , qµ)-channel were used for the flipped scenario, where the role of q is played by the strange quark. The extrapolations of the limits for the single-and pair production cross sections are compared to the corresponding theoretical resonant cross sections for V 1 in Fig. 4. The comparison for the case of V 3 is given in Fig. 5. We find good agreement with the similar extrapolations for the case of single-production in Ref. [79] . 3 We provide a list of possible mass reaches at future colliders for the leptoquark V 1 , in each of the three flavor scenarios, in Tab. III, for both pair-and single production channels. The reach for V 3 is given separately in parentheses, if different from the reach in V 1 .
As can be seen from Figs. 4 and 5 the theoretical predictions for single production in association with muons at 14 TeV and 27 TeV colliders turn out to be rather small, below the projected sensitivity. In Fig. 8 we compare the expectations for the single production cross sections for the flavor scenarios and different future collider experiments; the values span up to two orders of magnitude.
Observation of a single production signal with a cross section that is much larger than those shown in Fig. 8 would point to a leptoquark that is unrelated to R K,K * , since the simultaneous leptoquark couplings to all three flavors (d, s, b) are restricted by the low energy flavor-changing-neutral current (FCNC) observables, such as kaon decays.
In Fig. 7 we compare the cross sections for single-and pair production in the (M V , √ s) plane, where √ s denotes the center-of-mass energy of the p-p collisions. Pair production is instrumental for the discovery or the exclusion of vector leptoquarks in the region of a few TeV. For large masses and scattering energies, the cross sections of the single production turn out larger than those of the pair production -the corresponding regions are located on the right of the solid red lines in the first three plots of the Fig. 8. Notice that the cross sections of the single production vary significantly with the different flavor scenarios. In case of a signal discovery, b-tagging would be important in order to confirm the connection to the R K,K * -anomalies.  (14). In the regions to the right of the red lines the single leptoquark production cross section is larger than the pair production cross section. In the regions to the right of the blue lines (up to MV 1 ∼ 10 TeV) the resonant leptoquark production cross section is larger than the single production one, see text. In the plot to the lower right we show in addition the pair production cross section σ(pp → V +2/3 1 V −2/3 1 ). All plots are for κ = 0.

IV. CONCLUSIONS
Leptoquarks are flavorful -a feature that allows for rich phenomenology and model-dependence alike. While the recent evidence reported by the LHCb Collaboration [4] for the breakdown of lepton universality in rare semileptonic b-decays has yet to be confirmed by experiments and in other observables, taking the data at face value provides informative directions in the leptoquarks' parameter space: Among the spin 1 leptoquark representations only V 1 and V 3 induce sufficiently large contributions explaining R K,K * at tree-level [14], with coupling over mass ratio fixed (1). Here we study the vector leptoquark reach at the LHC and beyond targeting this parameter space. Specifically, we analyze signatures from leptoquarks with couplings to second and third generation quark doublets, and to muons. The reason why this simplified framework is sensible is two-fold: There is presently no necessity to consider couplings to electrons, and flavor symmetries explaining neutrino masses and mixing result in leptoquark couplings to a single lepton species [21]. We stress that dedicated searches for leptoquarks decaying to leptons other than muons [23] are well-motivated and complementary, however, beyond the scope of this work.
We work out single-and pair production cross sections in three quark flavor benchmarks: a hierarchical one (10), with dominant coupling to third generation quarks, a flipped one (12), with dominant coupling to second generation quarks, and a democratic one (13). Reinterpreting a recent ATLAS search for pair-produced scalar leptoquarks [22] we obtain the mass limit M V1 > 1.4 TeV for κ = −0.3, and higher otherwise. Analogous limits can be derived for V 3 , neglecting mass splitting within the multiplet, as M V3 > 1.6 TeV for κ = −0.3. Limits for gauge-type leptoquarks (κ = 1), or without the κ-term are stronger, see Sec. III B.
The  Table III. For κ = 1, the maximal reach for V 1 is 3 TeV (HL-LHC), 5.5 TeV (HE-LHC) and 19.9 TeV (FCC-hh). The reach for the triplet V 3 is similar in single production and the hierarchical pattern, and generically larger otherwise. All cross sections become larger for larger value of the parameter |κ|, as illustrated in Fig. 2, and improve the mass reach. Results are based on extrapolations of CMS [25] and ATLAS [22] searches. Pair production has larger cross sections due to the strong interaction until phase space suppression kicks in and single production takes over, as demonstrated quantitatively in Fig. 7. Single production is, however, valuable on its own as it is sensitive to the flavor patterns in the new physics sector. The flipped pattern with subject to the larger pdf gives largest cross sections, followed by the hierarchical one, see also Fig. 8. Recent works suggest to study resonant production from lepton-quark fusion at the LHC via lepton pdfs. Similar to single production the cross section is sensitive to flavor, as shown in Fig. 6 for the HL-LHC. A full study of efficiencies also at future machines is beyond the scope of this work.
Note also that the patterns (10), (12), (13) are simplified and in general lepton flavor violating signatures can arise in leptoquark decays, e.g., [14]. Allowing for significant entries ( * ) in the patterns would open up further search channels, and reduces leptoquark branching ratios in the signal channels studied here. Note that this rescaling effect is linear in single production, and quadratic in pair production, but leaves the qualitative features, such as quark flavor hierarchies of our analysis intact. A study in concrete, full flavor models [13] is beyond the scope of this work.
We conclude that leptoquark searches at the LHC are very well motivated by flavor physics, although covering the full mass range supported presently by rare processes requires higher energies. The actual confirmation of the R K,K *anomalies would strengthen the case for a corresponding machine. Observing leptoquarks directly would disentangle mass from couplings, as in (1), and could distinguish flavor patterns. Besides being striking signals from beyond the standard model, this would allow to make progress towards the flavor puzzle.
In order to estimate the cross section of singly produced V 3 with signature pp → jµµ in the flipped and democratic scenarios, we use σ(pp → V where the factor of 2 stem from adding the CP-conjugate of the process. Assuming that the pdfs for the strange and charm quarks are roughly the same, one obtains σ(sg → V 2/3 3 where the ( √ 2) 2 is an isospin factor (4). Using the different branching ratios given in Tables I and II,  (pp → jµµ), whereas an explicit evaluation results the ratio to be 2.5, and dropping for larger masses to 1.5, for relevant ranges of ( √ s, M V ). Numerical results from the explicit evaluations of σ(cg → V