Lepton-Quark Fusion at Hadron Colliders, precisely

When a TeV-scale leptoquark has a sizeable Yukawa coupling, its dominant production mechanism at hadron colliders is the partonic-level lepton-quark fusion. Even though the parton distribution functions for leptons inside the proton are minuscule, they get compensated by the resonant enhancement. We present the first computation of higher order radiative corrections to the resonant leptoquark production cross section at the Large Hadron Collider (LHC). Next-to-leading (NLO) QCD and QED corrections are similar in size but come with the opposite sign. We compute NLO $K$-factors for a wide range of scalar leptoquark masses, as well as, all possible combinations of quark and lepton flavors and leptoquark charges. Theoretical uncertainties due to the renormalisation and factorisation scale variations and the limited knowledge of parton distribution functions are quantified. We finally discuss how to disentangle the flavor structure of leptoquark interactions by exploiting the interplay between different production channels.


Introduction
Leptoquarks (LQs) are hypothetical new bosons that convert quarks into leptons and vice versa. The discovery of a leptoquark would represent a major breakthrough in our understanding of particle interactions, pointing towards an underlying quark-lepton unification at short distances. The phenomenology of TeV-scale leptoquarks is a very rich and mature subject, for a recent review see Ref. [1]. Leptoquarks at the TeV-scale are consistent with the non-observation of proton decay and can be found in wildly different settings beyond the Standard Model (SM). For example, they are in the spectrum of low-scale quark-lepton unification models à la Pati-Salam (see e.g. [2][3][4][5][6][7][8][9][10]). TeVscale leptoquarks also appear as pseudo-Nambu-Goldstone bosons of a new strongly interacting dynamics possibly related to the origin of the electroweak symmetry breaking (see e.g. [11][12][13][14][15]), or as a consequence of R-parity violation in supersymmetry (see e.g. [16][17][18]). On the one hand, they lead to distinct indirect modifications of low-energy flavor transitions, neutrino properties, top quark, electroweak precision, and Higgs physics. On the other hand, the direct production of a leptoquark at the LHC leaves a remarkable signature in the detector. Namely, a leptoquark would appear as a resonance in the invariant mass of a lepton and a quark jet.
Leptoquarks are colored just like quarks. Therefore, they are copiously produced in pairs in proton-proton collisions at the LHC by strong force [19][20][21][22][23][24][25]. A representative Feynman diagram is shown in Fig. 1 (a). In the limit of a small leptoquark coupling to quark and lepton (y q ), the scalar leptoquark production at hadron colliders is determined entirely by the strong coupling α s and the leptoquark mass m LQ . The phenomenology becomes more interesting once y q is increased. This is particularly relevant when establishing a connection with the low-energy flavor physics. The present indirect constraints on a TeV-scale leptoquark suggest that y q flavor matrix has a peculiar structure with some entries left unconstrained, and therefore possibly large. Taking a different perspective on the current data, in order to explain the existing experimental anomalies in B-meson decays [26][27][28][29][30][31][32][33] or muon g − 2 [34], some leptoquark couplings are required to be large.
If leptoquarks are indeed behind the origin of these discrepancies, there will be other production mechanisms beyond the QCD-induced pair production. To begin with, for a sizeable y q , there is an additional contribution with t-channel lepton exchange in qq fusion [35,36]. However, the production of two leptoquarks becomes quickly phase-space suppressed with increasing leptoquark mass. Therefore, often discussed in the literature is the single leptoquark plus lepton production from quark-gluon scattering [23,[37][38][39]. A representative Feynman diagram is shown in Fig. 1 (b). The production cross section for this process is proportional to |y q | 2 , but suffers less phase-space suppression. For a heavier leptoquark and a larger coupling, this production mechanism starts to dominate over the pair production. In this work, we are interested in a sizeable (yet perturbative) coupling range (i.e. 0.1 y lq √ 4π depending on the quark flavor), for which the production of a single leptoquark plus lepton becomes comparable or even favorable.
For example, Fig. 2 shows the relative comparison of different channels in the mass versus coupling plane when the leptoquark couples to down quark (left panel) and bottom quark (right panel). The upper edge of the vertical axis is chosen such that the t-channel induced pair production is suppressed compared to the pure QCD contribution. Nonetheless, the single leptoquark production plus the charge-conjugated (c.c.) process, dominates over the pair production in the large portion of the parameter space shown in Fig. 2. Relevant information on these parameters can also be extracted from indirect leptoquark effects at high-p T , such as Drell-Yan tails [40][41][42][43][44][45][46]. These probe complementary parameter space compared to both single and pair production (see Section 4 in Ref. [23]). The collider phenomenology of TeV-scale leptoquarks had a new twist recently. The precise extraction of lepton parton distribution functions (PDFs) [47] based on the LUX method [48,49] (see also [50,51]) facilitated another leptoquark production mechanism, the resonant leptoquark production [52][53][54]. The tree-level Feynman diagram is shown in Fig. 1 (c). The production cross section for the direct lepton-quark fusion is also proportional to |y q | 2 , but suffers even less phasespace suppression than the single leptoquark plus lepton channel. The difference between the two is the absence (presence) of a high-p T lepton. Therefore, the resonant channel cross section is always larger as shown in Fig. 2. Interestingly, this applies to all combinations of quarks and leptons involved. The ATLAS and CMS collaborations have extensively searched for leptoquarks in pair production and a single leptoquark plus lepton channel [55][56][57][58][59][60][61][62][63][64][65][66][67], however, the resonant production was not considered so far. Nonetheless, the phenomenological collider simulation in Ref. [54] shows that the resonant channel has a potential to probe the uncharted territory of interest in the mass versus coupling plane. It is therefore of utmost importance for leptoquark hunters at the LHC to place the resonant production mechanism at the top of their to-do list.
In this paper we fill in the gap on the theory side. Leptoquark toolbox for precision collider studies [23] includes leptoquark pair and single production at NLO in QCD. The scope of this work is to precisely calculate the resonant leptoquark production cross at the LHC including for the first time higher order radiative corrections and quantify the uncertainties from the missing orders and limited knowledge of parton distribution functions. The main result of this paper are the resonant leptoquark production cross sections at the LHC at NLO QCD plus QED with the corresponding uncertainties. These are reported in Tables 1 and 2, together with the complete set of NLO K-factors reported in Figs. 8, 9 and 10. Interestingly, we find that NLO QED corrections are as important as QCD corrections. Along the way, we discuss the interplay between different production mechanism and propose methods to determine the quark flavor inside the proton from which the leptoquark was created. The present study is limited to scalar leptoquarks and will be extended to include vectors in the future. Radiative corrections in the scalar leptoquark models are not sensitive to the details of the ultraviolet completion, in contrast to the vector case [68][69][70].
The paper is organised as follows. In Section 2 we set up the framework and present com- Comparison of cross sections for three leptoquark production mechanisms at the LHC ( √ s = 13 TeV). Shaded regions show the parameter space in the leptoquark mass versus coupling plane for which the corresponding cross sections are > 0.1 fb. The pair production pp → S † QLQ S QLQ cross section is shown in black, while the single (pp → S QLQ )+ c.c. and the resonant (pp → S QLQ )+ c.c. cross sections are shown in blue and red, respectively. In the left (right) panel, the leptoquark interacts primarily with the down (bottom) quark. The lepton flavors in the resonant production are shown with solid (τ ), dashed (µ) and dotted (e) lines. The electric charge of S QLQ is set to Q LQ = 2/3, however, the difference is negligible for Q = 4/3. For consistency, all cross sections are computed at NLO QCD (plus NLO QED for the resonant process) with the same central PDF set LUXlep-NNPDF31_nlo_as_0118_luxqed (v2) [47]. The first two processes are computed using the leptoquark toolbox [23], while the resonant production is taken from Section 3. The pair production from the t-channel leptoquark exchange is negligible in this coupling range. pact analytic expressions for the relevant partonic cross sections stemming from loop calculations detailed in Appendices A and B. In Section 3 we perform a numerical calculation of the hadronic cross section for the resonant leptoquark production at the LHC using the most recent lepton parton distribution functions. Supplemental numerical results are left for Appendix C. We finally conclude in Section 4.

Scalar leptoquark resonant production
The inevitable condition for a field coupling quarks and leptons at the tree level is to transform in the (anti)fundamental representation of the SU (3) part of the SM gauge group. The interaction between leptoquark and gluons is then completely specified and forms the basis for NLO QCD calculations. In contrast, the electroweak part of the SM allows for leptoquark representations involving different SU (2) L × U (1) Y multiplets with the corresponding hypercharges Y . As we are interested in evaluating the NLO QED corrections, the only relevant information is that after electroweak symmetry breaking, the possible absolute electric charges for any component of the SU (2) L multiplet are |Q LQ | = {1/3, 2/3, 4/3, 5/3}, in the units of the positron charge. Therefore, to assess the NLO QCD plus QED corrections to the resonant leptoquark production, we can treat one component inside the multiplet at a time, with the production of other components corresponding to separate processes.
The fermion content is the SM one, and the quark-lepton interaction with the scalar leptoquark S Q LQ of charge Q LQ is given by where y L,R q are 3×3 matrices in flavor space, encoding the most general form of Yukawa couplings. The chiral fermionic fields q L,R and L,R (note the left-and right-handed chiral projectors P L,R ) correspond to charge and mass eigenstates after the electroweak symmetry breaking. (Fermion mixings when going from the interaction to the mass basis, both in the quark and lepton sectors, are already absorbed in Yukawa matrices y L,R q .) Depending on Q LQ , some fermionic fields in Eq. (2.1) are charge-conjugated from the usual SM definitions, for example L ⊃ −y R ue u L e C L S † 1/3 . We can use the proton composition to precipitate lepton-quark fusion involving quark flavors u, d, s, c, b and charged leptons e, µ, τ . When calculating partonic cross sections, we will work in the limit of disregarding all fermion masses, which is an excellent approximation given the energy of the collisions. Neutrinos are not created in photon splitting and cannot be generated inside the proton at the order we are interested in. In the absence of fermion masses, possible interference terms involving left-and right-handed Yukawa couplings vanish. This allows us to independently treat processes in which leptoquark is resonantly produced by the same flavor combination of quarks and leptons, but of the opposite chirality. Additionally, the resonant leptoquark production is specified by one entry in the chiral Yukawa matrix irrespective of all other entries. When several flavor couplings contribute to the production of the same leptoquark, the individual contributions to production cross section factorise, and we add them separately.
In full generality, we summarize that scalar leptoquarks (SLQs) are SU (3) triplets, with four possible values of the electric charge, and their resonant production cross section is determined by the entries in the Yukawa matrices without interference. This exhausts all possibilities for SLQs and we conclude that our computation could be easily matched to any model containing these particles. Moreover, we note that neglecting fermion masses causes all one-loop corrections proportional to Yukawa couplings to vanish. Accordingly, for the case of SLQ, the dominant NLO effects originate from QED and QCD. 1 The relevant NLO QED (QCD) corrections to partonic cross section are calculated in Appendix A (B). The hadronic cross section is obtained after convoluting the relevant partonic cross sections,σ, with the parton distribution functions, f i and f j , in the following way, where ξ = m 2 LQ /s, √ s is the collider center of mass energy, y is the fraction of proton momentum carried by the parton labeled by i, and z = m 2 LQ /ŝ, withŝ being the partonic-level center of mass energy. The sum goes over ij = {q , g , qγ}, with the individual cases corresponding to Eqs. (2.3), (2.4), and (2.6), respectively. 1 The situation is different in the case of vector leptoquarks (VLQs). The calculation of NLO corrections for these particles necessarily involves details depending on the UV completion that embeds them. For instance, in many popular extensions of the SM, VLQs are accompanied by the massive color octet affecting NLO QCD contributions to processes involving VLQ in a nontrivial way [68,69]. For the moment, we focus on the resonant production of the SLQs, while postponing a detailed analysis of spin-1 leptoquarks for future work.

Next-to-leading order QCD corrections
The hadronic cross section for resonant leptoquark production is set by the size of the colliding parton densities, and the size of the parton level cross section. The Yukawa couplings are O(1), and at the leading order (LO), the partonic cross section scales asσ 0 ∝ |y ql | 2 . The parton density for gluons and quarks can be viewed as a sum of terms n (α s L) n , where α s is the QCD coupling, L = log(µ 2 F /Λ 2 ), with µ F representing the factorisation scale, and Λ is the typical hadronic scale. The QCD coupling is evaluated at the factorisation scale and its size is set by α s ≈ 1/L. We conclude that gluon and quark PDFs are non-perturbative objects of O(1). In contrast, the photon density is the first order QED effect and its size is determined by αL n (α s L) n , where α is the QED coupling. Further, as a result of photon splitting, lepton PDFs are generated at the next order in QED and their size is given by α 2 L 2 n (α s L) n . We apply the same QED to QCD coupling comparison already employed in [47][48][49] and use that α ≈ α 2 s . Accordingly, in terms of α s , the size of the photon density is O(α s ), while the lepton densities are O(α 2 s ). The size of the LO hadronic cross section for resonant leptoquark production is then . Therefore, the typical QCD correction coming from O(α s ) diagrams represents the contribution to hadronic cross section which is O(α 3 s ). The virtual corrections from gluon loops ( Fig. 5) are summed with the diagrams involving the real gluon emission (Fig. 6) to obtain the IR safe partonic cross section where C F = 4/3 and µ F , µ R are the factorisation and renormalisation scales, respectively. The remaining O(α s ) diagrams that contribute to the resonant leptoquark production involve gluons in the initial state (Fig. 7). The partonic cross section in this case readŝ (2.4) where T R = 1/2. The MS scheme was utilized both for factorisation and renormalisation. The NLO QCD corrections are universal for all leptoquark types. More details about the partonic cross section calculation can be found in Appendices A and B.

Next-to-leading order QED corrections
The NLO QED corrections are provided by processes where the initial lepton is replaced by a photon splitting into lepton pairs (Fig. 4). We estimate the size of these corrections by α s power counting for the leptoquark production via γ + q → + LQ. When convoluted with the corresponding PDFs, the size of the resonant cross section is Interestingly enough, the QED corrections are of the same order as the typical QCD corrections and their inclusion is essential in assessing the NLO effects in resonant leptoquark production. Employing the MS factorisation scheme, the partonic cross section readŝ The logarithmic part is universal for all leptoquark types since it originates from photon splitting to charged lepton pair, while the charge dependence is encoded in the functions X Q LQ (z), where subscripts {1/3, 2/3, 4/3, 5/3} correspond to electric charge of the leptoquark. Since the CP is conserved, the same formulas hold for the charge-conjugated processes. Note that loop diagrams involving photons are higher order in QCD coupling. Namely, the size of these diagrams, for the 1-loop corrections to partonic-level cross section involving photons,σ q , is given by and we neglect them. The detailed derivation of the NLO QED corrections is presented in Appendix A.

Numerical results and discussion
We carry out a numerical calculation of the hadronic cross section for the resonant leptoquark production in pp collisions. We consider the most general flavor structure of the leptoquark coupling y q to a quark q ≡ d L , d R , u L , u R , s L , s R , c L , c R , b L or b R , and a lepton ≡ e L , e R , µ L , µ R , τ L , or τ R . All options for the leptoquark electric charge are considered, |Q LQ | = 1/3, 2/3, 4/3 and 5/3. Cross sections are calculated for every q combination separately assuming y q = 1. As a reminder, the total cross section is simply the sum over different channels, σ = q, |y q | 2 σ q . We compute the process and its charge conjugate at leading and next-to-leading order in QCD and QED. We scan over the large leptoquark mass window m LQ = [500 − 5000] GeV relevant for the future studies at the LHC. As a benchmark, the collider center of mass energy is set to √ s = 13 TeV. Partonic cross sections are convoluted with LUXlep-NNPDF31_nlo_as_0118_luxqed (v2) parton distribution functions derived in Ref. [47]. To this purpose, we employ the Mathematica package ManeParse [71] for manipulating the LHAPDF grids [72]. The PDF extrapolation in Q 2 is checked by solving the corresponding DGLAP equations using Hoppet [73] in accordance with the prescription from [47]. Also, the running of the gauge couplings with the renormalisation scale is appropriately included. The central renormalisation and factorisation scales are set to Table 1.
Inclusive cross sections in pb for the resonant leptoquark production from up-type quarks, pp → LQ + charge-conjugated process, as a function of the leptoquark mass m LQ at √ s = 13 TeV. The cross section σ S 1/3 (σ S 5/3 ) corresponds to the resonant production of scalar LQ with absolute electric charge 1/3 (5/3) when the associated Yukawa coupling strength is set to one, y q = 1. The second column denotes which quark-lepton pair couples to the corresponding leptoquark. First (second) uncertainty is due to the renormalisation and factorisation scale variations (PDF replicas), and is given in per cent units. See Section 3 for details.
We estimate the uncertainty from the missing higher order corrections by varying the scales in the range Independently of µ R and µ F scale variations, the renormalisation group running of the leptoquark coupling y q (µ) in the range µ ∈ [0.5 − 2] m LQ leads to the cross section prediction uncertainty of about 4% across the entire m LQ window. The uncertainties due to the parton distribution functions are calculated by the method of replicas [51,74]. In particular, we report the standard deviation of the result calculated over one hundred replicas as the PDF error.
The lepton and antilepton PDFs are numerically the same. However, this is not the case for the light quarks, implying that e.g. ue − induced cross section is different fromūe + . We therefore report the cross sections for pp → LQ + the charge-conjugated process in Table 1 Table 2. Inclusive cross sections in pb for the resonant leptoquark production from down-type quarks, pp → LQ + charge-conjugated process, as a function of the leptoquark mass m LQ at √ s = 13 TeV. The cross section σ S 2/3 (σ S 2/3 ) corresponds to the resonant production of scalar LQ with absolute electric charge 2/3 (4/3) when the associated Yukawa coupling strength is set to one, y q = 1. The second column denotes which quark-lepton pair couples to the corresponding leptoquark. First (second) uncertainty is due to the renormalisation and factorisation scale variations (PDF replicas), and is given in per cent units. See Section 3 for details. and Table 2 (down-type quarks), at NLO QCD + QED accuracy. Thanks to the inclusion of radiative corrections computed in Section 2, the uncertainties due to the {µ R , µ F } scale variations are at the level of few per cent for all leptoquark charges, as well as, quark and lepton flavors and benchmark masses. The uncertainties due to the parton distribution functions strongly depend on the quark flavor and the leptoquark mass. In particular, the total uncertainty becomes dominated by the limited knowledge of the heavy quark PDFs when m LQ is several TeV. Next-to-leading order K-factors, defined as the ratio of NLO to LO results, are shown in Appendix C in Figs. 8 (electron), 9 (muon), and 10 (tau) for all possible quark and lepton flavors and leptoquark charges. One notable example is shown in Fig. 3 in the main text. These are calculated using the central PDF set and the central scales µ R = µ F = m LQ for the LO cross section, while at NLO, we consider {µ R , µ F } scale variation with the central PDF set. The red (orange) bands are with (without) NLO QED corrections. 2 In all cases considered, the error band dramatically shrinks, illustrating the importance of the NLO QED corrections. Interestingly, both QCD and QED corrections are large, however, they partially cancel in the total cross section. Inspecting Figs. 8, 9, and 10 we conclude that K-factors typically exhibit only a slight dependence on the leptoquark mass and electric charge, as well as, lepton flavors. In this calculation, we sum up cross sections for the process pp → LQ and the charge-conjugated process before taking the ratio. We checked that the individual K-factors for the two are very close to each other, thus we report only the K-factors for the sum.
We also study the dependence of the NLO K-factors on the PDF uncertainties. In particular, for every PDF replica we compute σ NLO /σ LO . We then derive the 68% confidence level range around the central PDF prediction. Interestingly, this band does not exceed the NLO QCD + QED scale variation band, except for a very heavy leptoquark close to the edge of the considered mass range, where the PDF errors are O(1) for some flavors. In other words, PDF uncertainties cancel in the ratio to a good approximation. We therefore conclude that the K-factors reported in Appendix C are robust, and will not change significantly by more precise PDFs in the future.
On the practical side, the existing leading order generators are missing the leptonic shower crucial to properly simulate the resonant leptoquark events. However, this shortcoming will soon be resolved, see the third footnote in Ref. [54]. Once this is in place, the K-factors derived in this paper can be directly applied to the future LHC resonant leptoquark searches to correct the overall signal yield. The main experimental difference between the single leptoquark plus lepton production and the resonant leptoquark production is the p T spectrum of the accompanied lepton. In particular, the lepton is hard (soft) in the former (latter) case. Therefore, measuring the lepton (or the leptoquark) p T distribution will enable efficient discrimination between different leptoquark production mechanisms at the LHC. To this purpose, it is crucial to have a good theoretical control over the p T spectrum. Our Appendix could serve as a starting point for this calculation. Note, that the leptoquark searches so far required the presence of two charged leptons which effectively vetos the resonant mechanism.
The leptoquark signature is quite unique; it will show up as a resonance in the jet-lepton invariant mass distribution. To study the flavor structure of the underlying interactions, one can make use of the flavor tagging of the decay products. Unfortunately, on the quark side there is a big degeneracy among light quarks u, d, s and c which are somewhat distinguished from the b quark. The task is even more difficult on the production side. The ratio of the rates for the single leptoquark plus lepton production and the resonant leptoquark production does not depend on the value of the leptoquark Yukawa coupling, however, it is sensitive to the initial quark flavor. This can be used to determine the flavor structure of the dominant leptoquark coupling in production. We have checked that the ratio drops quickly with the leptoquark mass and the discrepancy is more pronounced for sea quarks than for valence quarks.
Another observable relevant for the leptoquark flavor physics in high-p T collisions is the ratio of the resonant rate pp → LQ to its charge-conjugated process. For heavy c and b quarks the two rates are the same, while for the valence quarks the two rates can differ by a factor of O(10). We have checked that this observable indeed has a discriminating power, however, a dedicated analysis is needed to make a quantitative statement. We also noticed a large PDF uncertainties in the prediction of this ratio attributed to the poor knowledge of sea quarks at large x. Therefore, the success of this method depends on the improvements in measuring sea quark parton distribution functions.

Conclusions
A discovery of a leptoquark at the Large Hadron Collider would fundamentally change our understanding of particle physics, pointing towards a microscopic theory where quarks and leptons unify. Viable extensions of the Standard Model with TeV-scale leptoquarks exist, and are safe on proton decay and dangerous flavor changing neutral currents. Moreover, these models have recently received a large attention within the community. Namely, leptoquarks in the TeV mass range provide an elegant explanation of the long-standing hints on the lepton flavor universality violation in B-meson decays, as well as, the anomalous magnetic moment of the muon.
Leptoquark collider searches so far were mainly focused on the pair production mechanism driven by QCD interactions, while the role of the defining leptoquark interaction to a quark and a lepton was invoked in decays. However, interesting flavor effects occur when the leptoquark coupling is large(ish) [75][76][77][78][79][80][81][82][83][84][85][86][87][88], consequently predicting richer collider phenomenology on the production side. Building on Refs. [47,54], in this paper we study the resonant leptoquark production mechanism. Namely, the quantum fluctuations allow for a small presence of a lepton inside the proton which fuses with a quark from the other proton, to produce a leptoquark. The smallness of the lepton distribution is overcome by the resonant enhancement, providing this mechanism with the largest cross sections of all when m LQ 1 TeV and y q ∼ O(1), see Fig 2. We calculate for the first time next-to-leading order QCD and QED corrections to the resonant leptoquark production at hadron colliders. The present study is limited to scalar leptoquarks while the vector leptoquark case is left for the future work. The total cross section is given in closed form in Eqs. (2.2), (2.3), (2.4), and (2.6), and the detailed derivation is carried out in Appendices A and B. This formula is numerically integrated with the most recent lepton PDFs [47] to obtain the hadronic cross sections at the LHC. The main numerical results are reported in Tables 1 and 2, and in Figs. 8, 9 and 10. The calculation is performed for a set of benchmark points in the mass range relevant for the future searches, as well as, for all possible lepton and quark flavors and leptoquark charges. Importantly, our results are applicable for a general scalar leptoquark model with arbitrary flavor couplings.
We find that both QCD and QED corrections are large and are of similar size. However, they come with the opposite sign and cancel out in the final cross section, leading to somewhat smaller corrections of the tree-level result than initially expected. However, the advantage of our calculation is that we are now in position to reliably estimate the theoretical uncertainties. On this note, we observed a dramatic reduction of the renormalisation and factorisation scale variation uncertainties after the inclusion of QED corrections on top of the QCD ones. This is nicely illustrated in Fig. 3 with the red band. The leading source of theoretical error at this point is the limited knowledge of the parton distribution functions, in particular, the sea quark PDFs at large x. The breakdown of different uncertainties is summarised in the predictions for the total cross sections in Tables 1 and 2. The complete set of NLO K-factors is reported in Appendix C and can be straightforwardly applied in the future experimental searches at the LHC for the most general leptoquark model.
Finally, should a leptoquark be discovered at the LHC, precision measurements of the resonant process and its charge-conjugate, as well as, the single leptoquark plus lepton production, would help to deduce the flavor of the leptoquark interactions. Hopefully, synchronised deviations would show up in the low-energy flavor transitions to confirm this picture.

Acknowledgments
We thank Gino Isidori, Javier Fuentes-Martín, Matthias König and Pier Francesco Monni for useful discussions. The work of AG has received funding from the Swiss National Science Foundation (SNF) through the Eccellenza Professorial Fellowship "Flavor Physics at the High Energy Frontier" project number 186866. The work of AG and NS is supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, grant agreement 833280 (FLAY).

A NLO QED corrections to resonant production
The QED corrections to the resonant leptoquark production correspond to processes involving a photon in the initial state splitting into a lepton pair. As explained in Section 2.2 , the inclusion of these corrections is necessary for calculating the resonant production at O(α 3 s ), which is a typical size of the NLO QCD corrections. Due to different electric charges, the QED corrections to production cross section will differ for various leptoquark types. We can resonantly produce all types of scalar leptoquarks using the suitable partons inside the colliding protons. The possible combinations areū along with the corresponding charge conjugated processes, where the fraction in the subscript denotes the leptoquark electric charge. The QED correction for each of the listed process is given by three diagrams shown in Fig. 4. The initial state quark and photon create the complementary scalar leptoquark together with a soft charged lepton in the final state.
Generically, the amplitude for the process γ(p 1 )+q(p 2 ) → l − (k)+LQ + (q) obtained by interfering diagrams reads where p 1 (p 2 ) is the four-momentum of the photon ((anti-)quark) in the initial state, while q (k) is the four-momentum of the leptoquark (soft charged lepton) in the final state. The fermionic wave-function ζ(p 2 ) could either stand for the particle or anti-particle, depending on the produced leptoquark type. The partonic cross section calculation presented is the same for all types of leptoquarks, with the only difference provided by different particle charges, Q q denoting the (anti-)quark and Q LQ the leptoquark charge. To express the kinematics, it is convenient to use the center of mass frame in which with z = m 2 LQ /ŝ. Additionally, the relation between partonic Mandelstam variablesŝ = (p 1 +p 2 ) 2 andt = (p 1 −k) 2 becomest = −ŝw(1−z), where w = (1−cos θ)/2. The collinear divergences that appear when soft lepton is emitted parallel to the photon are regulated using dimensional regularisation with d = 4 − 2 . Averaging over the initial, and summing over the final polarisations and colors, the averaged squared matrix element can be written as where d − 2 in the denominator counts the polarisations of the massless gauge bosons in ddimensions, and M 2 div (M 2 fin ) denotes the part of the averaged squared matrix element that will produce the IR-divergent (IR-finite) contributions to the cross section after integration over the phase-space. In terms of w and z, they can be written as Moreover, the integration over the 2-body phase-space in d-dimensions with the corresponding flux factor, expressed in terms of w and z, can be performed as 1 16πŝ The integral over M 2 fin is IR safe for d = 4. The result for the finite contribution to the partonic cross section iŝ On the other hand, the phase-space integration over M 2 div induces IR-poles in the partonic cross section. Regulating the integral we obtain . (A.10) The IR-pole inσ div becomes explicit after w −1− is expanded around = 0 to give a distribution with the plus distribution defined such that The contribution to the partonic cross section containing a collinear divergence is then Note that the divergence is universal for all combinations listed in (A.1) as (Q q − Q LQ ) = −1, for all of them, and the corresponding coefficient is identified with the leading-order photon-to-charged lepton splitting function P ←γ (z) = z 2 + (1 − z) 2 . (A.14) In order to calculate the measured hadronic cross section for these processes, we need to convolute the partonic cross section with the corresponding quark and photon PDFs. In the procedure, the collinear singularity can be absorbed into the bare PDFs at a factorisation scale µ F . Here, we utilize the MS factorisation scheme by adding the counter term to the partonic cross section, consistent with the MS prescription used in extracting parton density functions presented in [47]. Combining relations (A.9), (A.13) and (A.15), we find that the result is finite, factorisation scale dependent, and can be written in the following form The non-universal pieces for different scalar leptoquarks, after the electric charges are replaced, read (1 − z)(7 + 37z) + Adding the virtual contributions to the tree-level amplitude, the NLO amplitude may be written as The chiral Yukawa couplings result in the vanishing leptoquark contribution to the massless fermion wave-function, with gluon providing the only contribution where L IR(UV) µ = log(µ 2 F(R) /m 2 LQ ), with µ F and µ R denoting the factorisation and renormalisation scales, respectively, and C F = 4/3. Similarly, the leptoquark two-point function Σ LQ (q 2 ), with q being the leptoquark four-momentum, receives no contribution from fermions, and the only effect is caused by the leptoquark coupling to gluons. We renormalise the leptoquark mass on-shell, with the wave-function correction defined as Taking the on-shell limit for the resonant leptoquark production q 2 = m 2 LQ , the correction becomes while the vertex correction, evaluated on-shell reads Combining the individual contributions listed above and integrating the averaged matrix element |A NLO | 2 over the leptoquark phase-space, we obtain the virtual correction to the partonic cross section for the leptoquark resonant production z).

B.2 Real QCD corrections
The calculation of the real QCD corrections closely follows the steps described in Appendix A for the QED case. We note that the results of this calculation already exist in the literature [89,90] (see also [91]), which we have checked and found complete agreement. With that in mind, here we present the corresponding results, and commit to these references for more details. The first process to consider is the one with the soft gluon in the final state, which can happen either by emission from the quark or the leptoquark as shown in diagrams in Fig. 6.
(a) (b) Figure 6. Diagrams for the process q + l → g + LQ contributing to the resonant leptoquark production at O(α s ).
The partonic cross section for this process is given bŷ As expected, the inclusion of the gluon radiation provides the IR divergences that exactly cancel the ones present in the virtual contribution (B.6). The combined result readŝ The second process which we need to take into account is the one with the soft quark in the final state, corresponding to the diagrams in Fig. 7.
(a) (b) Figure 7. Diagrams for the process g + l → q + LQ contributing to the resonant leptoquark production at O(α s ).
The partonic cross section for this process is given bŷ (B.9) where T R = 1/2 is the appropriate SU (3) color factor. Due to the massless quark which can be collinear to the gluon in diagram (a) of Fig. 7, this process needs the inclusion of the MS counterterm. This is exactly the same situation we already studied in the process γ + q → l − + LQ, and the universality of the log(µ F ) terms for the two processes in (B.9) and (A.16), up to color factors, becomes evident.