All possible first signals of gauge leptoquark in quark-lepton unification and beyond

We study possible current and future low-energy signals of the gauge leptoquark in quark-lepton SU(4) unification \`a la Pati-Salam. Taking fully into account the freedom in the generation mixing between quarks and leptons, we compile a catalogue of observables which currently form a border of the excluded part of the parameter space - hot candidates for first signals of new physics. We also determine the sensitivity needed in order to inspect a currently allowed part of the parameter space for several other measurements which are not included in this catalogue. We improve older similar works on this topic by taking into account more (and more recent) experimental measurements and by scanning the parameter space more densely. Furthermore, we study in a similar manner the SU(4) models with a small number of generations of extra leptons. We also discuss the minimal number of leptons needed in order to alleviate the contemporary discrepancies in the neutral-current B-meson decays.


I. INTRODUCTION
The main goal of this phenomenological study is to list all possible smoking gun signals of the Pati-Salam leptoquark.

A. Quark-lepton unification and gauge leptoquark
Quark-lepton unification (QLU)à la Pati and Salam [1,2] is an old idea motivated by the equal number of lepton and quark families and their similar electroweak behaviour. Technically, QLU is based on extending the QCD gauge factor SU (3) C to SU (4) C and accommodating the quarks and leptons in common 4-dimensional representations: The most characteristic prediction of QLU is the existence of a gauge leptoquark (LQ) U 1 transforming as (3, 1, + 2 /3) with respect to the Standard Model (SM) gauge group G SM = SU (3) C × SU (2) L × U (1) Y . The LQ has the following interactions with the fermions from Eq. (1): Here i ∈ {1, 2} is an SU (2) L index;d R ,û R ,ê R ,ν R denote the family triplets of the same-charge fermions in the mass basis, e.g.d R = (d R , s R , b R ), and similarlyq L and * gedeonova@ipnp.mff.cuni.cz † hudec@ipnp.mff.cuni.cz (corresponding author)l L are in the mass basis of their T 3 L = − 1 /2 components. The 3 × 3 flavour matrices V L , V R , V ′ R and the LQ mass m U1 are free parameters of the theory. QLU fixes the g 4 coupling at the scale of SU (4) C breaking and restricts V L , V R , V ′ R to unitary patterns, i.e.
The interactions of U 1 conserve baryon and lepton numbers but always introduce lepton flavour violation (LFV) and lepton flavour universality violation (LFUV) -see Appendix A. Hence, the gauge leptoquark is not restricted by proton stability nor by searches for neutrinoless double-beta decay, while extraordinarily high mass limits stem from flavour phenomenology: assuming V L = V R = 1 1, the experimental bound BR(K 0 L → e ± µ ∓ ) < 4.7 × 10 −12 [6] implies m U1 ≳ 2000 TeV. However, the gauge leptoquark has different phenomenology with different forms of V L,R .

B. Literature overview
Studies of the U 1 leptoquark have gained popularity in recent years as it has been identified as an excellent candidate to account for the neutral-current as well as charged-current B-meson anomalies (e.g. [7,8]). The benchmark setup for accommodation of the B anomalies as identified in Ref. [8] can be written as where ξ is a positive O(1) number. Clearly, such flavour and chirality pattern is incompatible with the conditions arXiv:2210.00347v2 [hep-ph] 13 Aug 2023 in Eqs. (3). For this reason, most of the current studies employ chiral vector-LQ models, based on more complicated gauge groups or on complete abandonment of the gauge nature of the U 1 field. Despite its inability to account for the discrepancies in the B-meson decays, the gauge leptoquark in the QLU framework is worth a detailed and dedicated study as it is a common feature of many specific models. Several top-down studies have already been published in the last decades.
In 1994, Valencia and Willenbrock [9] considered the cases where V L = V R are permutation matrices, i.e. where each lepton is coupled to a single quark, and studied various two-body meson and tau decays. They found that apart from K 0 L → eµ, the gauge LQ mass was for some mixing patterns limited from below to 250 TeV by R e/µ (π + → l + ν) or R e/µ (K + → l + ν), or by BR(B + → e + ν) to m U1 > 13 TeV. At around the same time, Kuznetsov and Mikheev [10] considered various (semi)leptonic K and π decays and the µ → e conversion on nuclei, and cast inequalities employing m U1 and elements of quark-lepton mixing matrices, virtually taking the full freedom in the quark-lepton mixing into account, but still tacitly assuming V L = V R . Apart from BR(K 0 L → eµ) and R e/µ (K + → l + ν), important bounds have been found to stem also from BR's of K 0 L → l + l − , K → πµe and from coherent µ → e conversion on titanium nuclei. Needless to say, both analyses [9] and [10] are outdated nowadays due to new experimental data.
Concerning more recent works, Ref. [11] considered K 0 L → eµ and B 0 → eτ for general forms of V L,R but did not confront the obtained limits with other measurements. In Ref. [12], which is the 2012 update of [10], also the B factory results on B and τ decays have been included and the general case V L ̸ = V R has been considered. A specific form of V L and V R has been found for which the stated LQ mass limit was as low as 38 TeV. However, as pointed out in Ref. [4], this finding is invalid because the authors forgot to include the predictions for the µ − e + final state when studying BR(B 0 → µ ± e ∓ ) and BR(B s → µ ± e ∓ ).
Finally, Smirnov [4] considered all kinematically allowed leptonic decays P 0 → l + l ′− for P 0 = K 0 L , B 0 , B s and took fully into account the freedom in the fermion mixing by performing a scan. The global lower limit stemming from these processes was found to be m U1 > 86 TeV (5) and the corresponding forms of V L and V R were given. We have verified the computations by completely recalculating Ref. [4].

C. Outline of our work
The main goal of this work is to identify all observables which currently determine the gauge LQ mass limit for some form of V L,R = (V L , V R ). These observables are excellent candidates for future New Physics (NP) signals since even a small improvement in the precision of their measurement shall explore a yet unexcluded part of the parameter space of the model. Hence, we call them possible first future signals of the gauge LQ.
Clearly, this is a more ambitious aim than just finding the global LQ mass limit which is the main result of Ref. [4].
In the analysis, we focus especially on the following: • We attempt to take into account all relevant observables in which the signal of the gauge LQ in the foreseeable future might be potentially found.
To this end, we employ the Python package flavio [13] which is capable of calculating predictions for hundreds of observables.
• More recent measurements are included.
• No ad-hoc assumptions are made on the form of V L,R . Keeping in mind that there is no physically meaningful measure on the parametric space, the setups which might be labeled as fine-tuned scenarios or small parts of the parameter space are not dismissed.
Section II describes the model in more detail. In Section III, the technicalities of the calculations are presented. In Section IV, we present the results and discuss the potential of various relevant forthcoming experiments. Then in Section V, we analyze in a similar manner the SU (4) C models extended by several generations of left-and/or right-handed leptons. We briefly conclude afterwards. In the three appendices we provide some additional details concerning the lepton flavour group in LQ models, the physics of the Z ′ boson, and the optimization of the scanning procedure, respectively.

II. MODEL DETAILS
The SU (4) C gauge symmetry can be realized in a minimal way within the gauge group [3,5]. This symmetry might be an intermediate stage of a left-right theory based on the Pati-Salam group G 422 = SU (4) C × SU (2) L × SU (2) R [1,2]. The spontaneous symmetry breaking (SSB) of G 421 proceeds in two steps as The generators T 1 C , . . . , T 8 C of the unbroken part of the SU (4) C symmetry form SU (3) C , while the weak hypercharge is given by Y = 2/3 T 15 C + R. During the first step of symmetry breaking, massive gauge leptoquark U 1 and massive Z ′ arise; the W and Z bosons acquire mass during the second step in a SM-like manner. For further details we refer to [3,5] or to Appendix B.
The fermion sector consists of 3 generations of the fields in Eq. (1). Independent quark and chargedlepton masses can be achieved by using both 1-and 15dimensional scalar representation of SU (4) C [2,3,5]. Concerning the neutrinos, in principle one can assume that they are of Dirac nature [3,[14][15][16][17]. In such a case, the tiny neutrino masses are obtained as a difference of two parameters of the order of the top-quark mass. To avoid such a fine-tuning, one might call for some form of a seesaw mechanism. The traditional type-I seesaw (studied recently in Ref. [18] in this context) would require that the SU (4) C breaking scale is so high that the gauge LQ would have no measurable low-energy phenomenology; hence, this case is not of interest for us. Nevertheless, unification of quarks and leptons is possible even at a low scale when employing the inverse seesaw [5,18] instead. This model has been recently studied in context of the B-anomalies addressed by the scalar leptoquarks by Faber et al. [19,20] and also by Fileviez P. et al. [21,22] with mutually conflicting conclusions.
In accordance with the inverse-seesaw model, we assume heavy ν R in this study. Nevertheless, as we shortly discuss in Section IV, the results would be essentially identical also in the Dirac-neutrino case.
The interplay among the flavour matrices V L , V R , V ′ R introduced in Eq. (2) and the weak interaction matrices V PMNS and V CKM is illustrated in Fig. 1. Adopting the standard (single-phase) parametrization of V CKM and V PMNS , no complex phases can be removed from V L or V R . By expanding the SU (2) L structure in Eq. (2), the interactions of the SM fermions with the gauge LQ can be rewritten as whereν is a neutrino flavour triplet in the weak interaction basis and P L,R = (1 ∓ γ 5 )/2 are chirality projectors.
As the particular form of the V ′ R matrix is inconsequential for all the considered low-energy processes, the relevant dimensions of the parameter space are given solely by 2 × 9 angles or phases of V L,R and by m U1 .
We do not take into account any other BSM field in the model. Especially, the scalar sector is neglected. Note that the free parameters of the full renormalizable model [5] (or [3]) indeed allow for the regime in which the gauge LQ signals dominate over those of the scalars. Notice also that the interactions of Z ′ are flavour-diagonal and hence, its effects in flavour physics are suppressed. If there is no intermediate stage in the G 421 → G SM symmetry breaking, m Z ′ is of the same order as m U1 ; in such a case, Z ′ can be also safely neglected. For more details, see Appendix B.

III. METHODS
In what follows, by parameter space we mean the 18dimensional set of forms of V L,R . A parameter point is an element of this set, parametrized by angles and phases λ L,ij and λ R,ij as described in Appendix C.
We have employed two different approaches to investigate a chosen parameter point. The simplified approach, adopted from Ref. [4] and detailed in Section III A, served as a primary stage providing basic but yet coherent insight into the parameter space. Within the concept, the identification of interesting parts of the parameter space was quite straightforward since it makes use of simple analytical formulae for observable predictions. The more robust approach, described in Section III B, is more comprehensive, but also much less intuitive since it is based on numerical packages which we have used mostly as a black-box tool.
The former approach served also as an important crosscheck which enabled us to find and correct an error in the flavio package. 1 Hence, even tough the presented results are based solely on the latter, we still find it worthy to present also the first approach below.
In Section III C, we describe how the analyzed parameter points have been chosen.

A. Simplified approach
This approach directly follows Ref. [4]. There are several aspects about this procedure worth mentioning: 1. The effects of the U 1 leptoquark are taken into account at the tree level.
2. Four-loop QCD running of the induced effective operators is taken into account [23]. For simplicity, the effective operators are defined at the 100 TeV scale, regardless of the considered LQ mass.
3. SM contributions to the considered processes are completely neglected in the calculation. To highlight this approximation, the corresponding predictions for branching ratios are labelled by BR V . The 1 There was a bug in the expression for the K 0 L,S → e ± µ ∓ amplitudes. measured BR's of the decays which have been already observed (i.e., K 0 L → ee, K 0 L → µµ, B 0 s → µµ) are taken as limits on BR V . Such a rough approximation is meaningful due to large relative theoretical uncertainties for the SM amplitudes.
4. Ref. [4] has taken into account the branching ratios of P → l ± l ′∓ decays where P = K 0 L , B 0 , B 0 s and ll ′ corresponds to various kinematically allowed combinations of leptons and antileptons. In our work, also the leptonic decays of K 0 S are considered. The limits on B 0 d,s → e ± µ ∓ are updated [24]. 5. No processes with neutrinos are analyzed; the study holds for both situations with light or heavy righthanded neutrinos.
6. The masses of electrons and muons in the final state are neglected, as well as the indirect CP violation in the neutral kaon mass eigenstates.
7. For given V L,R , the LQ mass limit is determined as the maximum of individual limits obtained from the considered observables. The decay responsible for the strongest limit is considered to be the candidate for the future first signal of the LQ for the investigated form of V L,R .
The branching ratio for a decay with light leptons only is calculated by the following formula: where the formfactors are f K = 155.72 MeV, f B 0 = 190.9 MeV, f B 0 s = 227.2 MeV and m P = m 2 P /(mq + m q ) with q and q standing for the index of the valence antiquark and quark of P , respectively. The gluonic corrections to the pseudoscalar quark currents amount to R V K = 3.47 and R V B = 2.1 [23]. The lepton-flavourdependent factor is a sum over two different helicity combinations where for weak eigenstates while for the CP eigenstates, Here + and − relate to K 0 L and K 0 S , respectively, and a stands for either a LR or a RL . See Fig. 2 for an illustration.
For processes with a single τ lepton in the final state, the expression for BR V in Eq. (9) must be multiplied by a phase space factor (1 − m 2 τ /m 2 P ) 2 . Along with that, the replacement (V L,R ) rτ → (V L,R ) rτ − (V R,L ) rτ m τ /(2m P R V P ) for r = q,q is applied in Eq. (11). For τ + τ − in the final state see Ref. [4].

B. More robust approach
In parallel with the previous approach, we have also performed a similar analysis using the family of generalpurpose open-source tools wilson [25,26], flavio [13,27], and smelli [28,29]. We present the features of this approach as a list which can be compared with that in the previous section.
1. The LQ interactions are matched onto the Standard Model effective field theory (SMEFT) at the tree level (similarly to the previous approach), yielding non-zero Wilson coefficients which multiply the following effective operators (with flavour indices suppressed): We have implemented a python function taking V L,R and m U1 as input arguments and returning a dictionary of SMEFT Wilson coefficients in the format compatible with the wcxf standard [30,31], which is used by the packages mentioned above.
2. The renormalization group (RG) running of the SMEFT effective operators from the scale µ = m U1 to the electroweak scale, the tree-level matching onto the Weak effective theory (WET) and further evolution to the meson-mass energy scales is handled automatically by the wilson package. The full numerical solution to the one-loop SMEFT RG equations (the 'integrate' option) is performed since we have exemplified that the 'leadinglog' approximation leads to O(1) relative differences in certain predictions. Analytical solution to the oneloop QCD and QED running equations is applied under the electroweak scale in wilson. For more details see [25] and references therein.
3. The SM contributions to the amplitudes of the calculated processes are automatically taken into account by flavio. As a result of this (and of the RG running), the predictions do not scale uniformly as m −4 U1 , which was a simplifying feature of the previous approach [see Eq. (9)].

The global likelihood tool smelli is employed.
This package uses flavio for predictions and confronts them with the measurements, including correlations. By default, version 2.2.0 of smelli takes into account hundreds of observables, most of which are, however, irrelevant for our scenarios.
On the other hand, the very interesting processes BR(B 0 d,s → e + e − ) as well as µ → e conversion on nuclei were not included. To this end, we have modified the smelli package to calculate also these observables.
The complete list of considered observables can be found in [32] or inferred from [33].

5.
No light right-handed neutrinos are assumed.

Light lepton masses are taken into account in
flavio for all observables, but indirect CP violation in neutral kaons remains neglected in the K 0 L,S → ll ′ decays. 7. For V L and V R fixed, we find m U1 for which the global log-likelihood calculated by smelli worsens by 4 units with respect to the SM. That value defines the lower LQ mass limit for this particular case. Then, the corresponding candidate for the future first signal of NP is the observable for which the individual pull between theory and experiment worsened the most compared to the SM case; the pulls have been obtained via the obstable method provided by the smelli package [28].
We have also tried different (more complicated) criteria, supposed to underpin scenarios in which the likelihood actually improves, but we ended up with qualitatively identical results.

C. Analyzing the parameter space
The final analysis has been performed within the more robust approach where analyzing a single parameter point typically takes over a minute on a usual computer. Apparently, an 18-dimensional parameter space cannot be rigorously explored just with a blind numerical scan. To this end, we have addressed the issue in two mutually complementary ways: • A series of random numerical scans has been performed, using a naïve measure Π ij dλ L,ij dλ R,ij , where i, j run only over the unfixed λ's. The gradual fixing of λ's proceeded along the following lines (the details can be found in Appendix C): Illustration of lower limits on the gauge LQ mass stemming from several observables, calculated along a onedimensional cut of the parameter space described in Appendix C. On this slice of the parameter space, the bounds are given by BR(K 0 L → e ± µ ∓ ) and BR(B 0 → e ± µ ∓ ). The mass limits are obtained using the approach described in Section III A.
In the first stage, 10 3 parameter points have been obtained with none of the λ's fixed. In majority of cases, the limiting processes were BR(K 0 L → eµ), BR(K 0 L → ee) and CR(µ → e, Au), i.e., the coherent conversion rate of µ → e on nuclei.
Following Ref. [4] in the second stage, about 2×10 3 parameter points have been obtained on the parameter subspace defined by BR V (K 0 L → ll ′ ) = 0, achieved by Eqs. (C10) and (C11). This is motivated by exploring the "steep valleys" on Fig. 3. Now CR(µ → e, Au) dominated almost all cases.
• We have compiled a list of relevant observables discussed in the recent review in Ref. [34] and investigated if they might become the future first signal.
For each of these observables, we have found either a parameter point for which this observable is the first future signal indeed, or an argument that such a point should not exist. A thorough effort has been made to include various special parts of the parameter space in the considerations.
Combining those two methods enables us to claim with a higher level of confidence that the catalogue in Table I is complete.
Examples of Feynman graphs underpinning the possible first signals of the U1 gauge leptoquark.
TABLE II. Examples of processes which are not listed in Table I. The third column shows predictions obtained during the numerical scanning following from the forms of VL, VR, and mU 1 which are fully compatible with all the current experimental limits. We also list the SM predictions for comparison.

IV. RESULTS
Tables I and III present the catalogue of observables which currently give the most stringent constraint on m U1 for some configuration of V L,R . These observables correspond to the future first signals as defined above. To fully appreciate the result, notice that even a very small improvement in precision of any experimental limit listed in Table I will probe a so-far allowed part of the parameter space of the model, and could potentially detect a NP signal -the only exception is the observed decay K 0 L → µµ for which the theoretical uncertainties within the SM dominate.
Conversely, under a very idealized assumption that the experimental sensitivity will grow uniformly for all the observables considered, no other observable could become the first observed signal of the gauge LQ. More realistically, the measurement precision of any other observable needs to be improved by a larger step in order to put a new constraint on the model parameters or to have a theoretical chance of observing a signal of the gauge LQ. How large these steps must be is shown for several important examples in Table II. A. Global mass limit -comparison with Ref. [4] As noted earlier, the simplified approach of Ref. [4] described in Section III A leads to the global lower leptoquark mass limit of 86 TeV. The corresponding V L,R is shown in the last line of Table III. However, when taking into account more observables in the more robust approach, m U1 = 86 TeV for this parameter point turns out to be in conflict with the bound on CR(µ → e, Au) by 3 orders of magnitude.
Nevertheless, we have found a form of V L,R which allows essentially the same mass (90 TeV, see Table III) even when all the constraints included in smelli are considered.

B. Possible first signals
Concerning searches for LFV, Table I contains limits on K 0 L , B 0 d,s → eµ and on the µ → e coherent conversion on nuclei; further searches for these processes are therefore of great interest. The remaining observables in Table I are all related to the leptonic decays of pseudoscalar mesons which are chirality suppressed in the SM and can be understood as tests of LFUV in the SM : Firstly, significant deviations could arise in the ratios of charged current decays R e/µ (P + → lν) with P = π, K when the LQ couples mostly to the electrons. Although the decay widths involved cannot be measured with the precision similar to the rare decays above, the deviations from the SM can be significant due to the interference among the NP and SM amplitudes. Subdominant contributions arise also from the other neutrino species as well as from the l L ν R final state if the right-handed neutrinos are light enough.
Secondly, limits on m U1 stem also from the observed BRs of K 0 L → ee, µµ and B 0 d,s → µµ. Concerning K 0 L → µµ, the experimental precision is better than the theoretical error estimates in the SM stemming from long-distance contributions [49,50].
Finally, a very interesting limit on the U 1 mass for some patterns of quark-lepton mixing is set by the recent LHCb search for K 0 S → µµ; the anticipated discovery of this decay after the upcoming LHC runs thus provides an exciting opportunity for the Pati-Salam-type leptoquark.

C. Other observables
In Table II, the P 0 → ll ′ decays that currently do not pose the most stringent bound on m U1 are listed, together with the predictions based on the parameters fully compatible with all the current experimental searches. All generated parameter points have been included.
As τ leptons are generally experimentally hard to handle, all processes involving τ 's belong to this category. In fact, 3 ∼ 4 orders of magnitude improvements in limits on B 0 d,s → lτ would be necessary in order to compete with the other constraints, which is far below the prospected sensitivity of Belle II [51] and hardly achievable even at LHCb at the high-luminosity phase. Furthermore, as explained in Appendix C, due to the unitarity of V L,R , the LQ amplitudes mediating of B 0 d,s → τ τ are severely limited by the probes of K 0 L → ll ′ and, thus, our predictions for the former essentially coincide with the SM. Hence, the expected sensitivity of Belle II at about 10 −6 for BR(B 0 → τ τ ) [52] shall not be an interesting probe of the considered model.
On the other hand, the experimental sensitivities to K 0 S , B 0 , and B s decays to e + e − require less than 1 order of magnitude improvement in order to probe the currently unexplored parts of the parameter space. Note [4]; currently, the muonic channel is measured more accurately. However, when further searches for NP in the P 0 → µ + µ − decays become limited by the SM uncertainties, new searches for B 0 d,s , K 0 S → ee will become essential.
No experimental limits on the decay K 0 S → eµ are available [38]. Comparing with the current limits on K 0 S → ee [44] and K 0 S → µµ [39], we reckon the required experimental sensitivity around 10 −10 for K 0 S → eµ might be reachable by KLOE II or LHCb.
Semileptonic decays like B → Kµµ or loop processes such as µ → eγ might become the dominant signals of chiral leptoquaks but not of the gauge LQ in the considered model as it inevitably introduces sizable Wilson coefficients C ℓedq which are experimentally more constrained (see, e.g., Ref. [53]). This part of our work is devoted to more complicated models featuring the vector leptoquark U 1 . Although they could be considered as aesthetically less appealing, such models have been studied thoroughly in the recent years, mainly due to the attempts to accommodate the B-meson anomalies. Generally, several tricks to circumvent the theoretical requirement of unitarity of V L and V R have been suggested in the literature. They can be divided into three categories, according to the paradigm abandoned: 1. Adding extra generations of fermions while maintaining the gauge symmetry group G 421 or G 422 [54][55][56][57].
2. Assuming more complicated gauge structure. Especially, the models based on the G 4N 21 = SU (4) C L × SU (N ) C R × SU (2) L × U (1) R gauge symmetry have become popular; here N = 3 or 4 and the QCD generators are given by In the basic setting of chiral quark-lepton symmetry [58,59], the left-handed fermions are charged by SU (4) C L while the right-handed ones transform non-trivially under SU (N ) C R . Hence, the U µ 1 field interacting with the left-handed quark-lepton currents is a chiral leptoquark -it has no or suppressed couplings to the right-handed currents, avoiding the scalar-type effective operators O ℓedq which are responsible for all the most stringent limits in Table I.
3. Assuming that the vector LQ is not a gauge field but a composite resonance formed by some more fundamental strongly interacting fields [68][69][70].
This work is focusing solely on the first option. Since the SM leptons do not entirely stem from the same SU (4) C representations as the quarks, we shall not use the term quark-lepton unification for these theories but rather call them extended SU (4) C models.

A. Specification of the models
Like in the previously considered SU (4) C × SU (2) L × U (1) R scenarios, see Eq. (1), the models contain 3 generations of each of the following chiral fermion SU (4) C quadruplets: Notice that we have slightly updated the notation by adding an "isotopic index" to the leptons living inside the quadruplets. On top of that, k L generations of SU (2) Ldoublet vector-like fermions and k R generations of weak-singlet vector-like fermions are assumed. Being SU (4) C singlets, these new fields are intact to interactions of the gauge LQ. After the G 421 → G SM symmetry breaking, they can mix with the leptons from the quadruplets. We assume that the 3 lightest eigenstates correspond to e, µ, and τ , while the k L + k R remaining ones are too heavy to be observed. As the weak hypercharges of the 3 known leptons are quite precisely measured, they must be composed solely from the fields 1 e R , 4 e R , 1 ℓ L and 4 ℓ L . For all practical purposes, it is sufficient to assume the following mixing pattern in the charged-lepton sector: where, generally, ℓ − denotes the electrically charged component of an ℓ doublet (notice that 1 e L ̸ = 1 ℓ L − and 1 e R ̸ = 1 ℓ R − ),ê =ê L +ê R is the triplet of light leptons while E R and 1 ℓ R − with their chiral counterparts E L and 1 e L form the heavy mass eigenstates. The form of the mixing in the heavy-lepton sector is irrelevant for our considerations. The blocks V e L and V e R are arbitrary unitary matrices of dimension 3+k L and 3+k R , respectively. Including the "non-standard" fields 1 ℓ R and 1 e L into the model ensures the ABJ anomaly cancellation and enables one to write down arbitrarily large Dirac mass terms for the vector-like pairs.
The Q = 0 components of 4 ℓ L and 1 ℓ L naturally follow their charged SU (2) L partners during the mixing at the first stage of SSB: those belonging to E L become equally heavy while the companions ofê L become the light neutrinos, eventually gaining mass after the electroweak symmetry breaking.
There are no extra quarks in the models and the transformation from gauge to mass eigenstates is given by 3×3 unitary matrices: Finally, let us have a look at the gauge LQ interac-tions. Like in previous sections, we assume that ν R are heavy due to the inverse seesaw [5,71] and therefore their interactions with the U 1 leptoquark are unimportant for the low-energy phenomenology. Interactions of U 1 with the other fermions can be rewritten as follows: The L L field on the last line is the heavy SU (2) L doublet containing E L as a component. Apparently, the novelty of such extended SU (4) C models consists in the fact that the unitary matrices V L,R , defined by the last line of Eq. (21), are now of dimension 3 + k L,R . Using the block-form notation only the 3×3 submatrices V 0 L,R are relevant for the interactions among the SM fermions. The larger the numbers k L,R of extra lepton generations, the more parametric freedom in V 0 L,R is available. With k L = k R = 3, one can already choose any form of g4 m U 1 V 0 L,R , which is all that is relevant for the low-energy phenomenology at the leading order, cf. Eq. (13). Similar models have already been studied in the literature, usually considering the cases equivalent to (k L , k R ) = (3, 0) [72], (0, 3) [55] or (3, 3) [54]. In this work, we focus on the more economical models with k L,R < 3, which are less challenging if one aims to capture all the possible NP signals in the model, but more restrictive if parameters leading to a chosen signal (such as the b → sµµ anomalies) are searched for.
Note that enlarging the dimension of V L,R is indeed the only practical consequence of extending the theory of QLU from previous sections: we assume that the extra leptons are too heavy to be observed and ignore the details of the scalar sector responsible for the mixing. A construction of the scalar sector leading to a chosen form of V L,R in similar models can be found, e.g., in Ref. [63].
Note that although we keep neglecting the Z ′ in the model, it may actually be relevant in some cases. The discussion of this issue is deferred to Appendix B.

B. First signals of gauge leptoquark in extended SU (4)C models
We have performed an analysis similar to that described in Section III for the extended SU (4) C models with (k L , k R ) = (1, 0), (0, 1), (2, 0), (0, 2), and (1, 1). Some details about the scanning procedure can be found in Appendix C.
With growing number of free parameters, more couplings can be "rotated away" from V 0 L,R to the other parts of V L,R . New interaction patterns become allowed, with lower lower limits on m U1 . Naturally, the catalogue of the first future signals (the observables which currently constrain m U1 for some form of V L,R ) grows with the growing dimensions of these unitary matrices. The results are captured in Table IV.
While a lot of effort has been spent to fully explore the parameter space in the cases (k L , k R ) = (1, 0) or (0, 1), the number of parameters for k L + k R = 2 is quite high and we admit that the corresponding lists in Table IV may not be complete.

C. Addressing neutral current B anomalies
During the last decade, several discrepancies in both charged-current and neutral-current B-meson decays have been reported [73][74][75][76][77]. Plenty New Physics interpretations have been suggested (see, e.g. [8,78]), including the U 1 leptoquark. Achieving the setup form Eq. (4) is meaningless within our restricted model as it requires so low scale of SU (4) C symmetry breaking that neglecting the other BSM fields would be inadequate. Nevertheless, reasonable considerations can be made once only the accommodation of the neutral-current anomalies is sought for. These anomalies include the tests of lepton flavour see Table I BR universality Table V). Further measurements indicate that the NP effect is in the b → sµµ channel [76,77].
In order to ascribe this effect to the gauge LQ, the elements V Lsµ and V Lbµ need to be non-negligible. To avoid the scalar-type operators O ℓedq involving electrons or muons, which are responsible for the most severe constraints found in Section IV, k R = 2 generations of extra leptonic SU (2) L -singlets are required; the model with dim (V L ) = 3 and dim (V R ) = 5 allows for the following setup: Note that a similar pattern for V L has been suggested in Ref. [55] and also in Ref. [59] within the SU (4) C L × SU (4) C R × SU (2) L × U (1) R framework where the couplings to the right-handed fermions are suppressed globally.
Adopting Eq. (24), the maximum likelihood fit is close to the simple case which improves the global log-likelihood function of smelli [28] by more than 14 units compared to the SM, i.e. log(L/L SM ) ≈ 14. Such a scenario accommodates well the R K ( * ) anomaly and also significantly mitigates the tension in the additional b → sµµ observables.
Using the standard normalization factor N = 16π 2 for the effective four-fermion operators in the weak effective theory at the 5 GeV scale, Eqs. (24) and (25) imply the following contributions of New Physics to the Wilson coefficients: C NP 9 bsee = C NP 9 bseµ = +0.24 , (27b) In comparison, the benchmark one-dimensional effective scenario with only C NP 9 bsµµ = −C NP 10 bsµµ = −0.53 [8] improves log-likelihood to log(L/L SM ) = 18; the simplified vector LQ setup in Eq. (4) leads to log(L/L SM ) = 30 as it also accommodates R D ( * ) . Note that the discussion in terms of confidence levels would be pointless since these models differ in number of free parameters.
Predictions for several important observables following from Eqs. (24) and (25) are given in Table V. As outlined in Section III B, the LQ has been integrated out at the tree level and the calculated LFV dipole operators responsible for µ → eγ arise solely from the one-loop RGE running of the Wilson coefficients. Thus, the predictions for the loop processes should be interpreted with caution.
In the scenarios with nonzero couplings V Lse , V Lbe , V Lsµ , and V Lbµ , the strongest bounds arise from B + → K + µ ± e ∓ and from the the LFV loop processes like µ → eγ (see Ref. [84] for a dedicated study). Generally, the constraints from the latter are quite strong. However, in the chiral leptoquark models with unitary interaction matrix, µ → eγ is suppressed by an analogue of the GIM mechanism. As the only non-vanishing element of V 0 R in (24) is essentially irrelevant for µ → eγ, the same applies also to our case. Note that Ref. [84] did not consider the subleading terms and hence found exactly zero contributions to µ → eγ for the case V Lse V L * bµ . Ref. [59] considered the case equivalent to V L from (24) and V R = 0, finding the constraint m U1 > 10 TeV based on the BaBar search [85] for B → Keµ. The very recent measurement by LHCb [82] has pushed this limit to 17 TeV for the considered interaction pattern.
Finally, let us note that although the Z ′ interactions are not lepton-flavour universal, the couplings in the particular case of Eq. (24) are lepton-flavour diagonal and, hence, the Z ′ does not mediate any flavour violating processes (see Appendix A for more details about lepton flavour). At the same time, with the mass around 20 TeV, Z ′ is also safely hidden to the high-energy searches at LHC. We elaborate on Z ′ in Appendix B.
To conclude, the interactions of the SU (4) C gauge leptoquark in a model with two extra weak-isosinglet charged leptons can accommodate the neutral-current Bmeson anomalies to a large extent. The suggested scenario can be excluded by future negative searches for B → Keµ at LHCb or Belle II.

VI. CONCLUSIONS
We have studied the phenomenology of the gauge leptoquark model with SU (4) C symmetry of the Pati-Salam type, taking into account the most recent experimental data. The catalogue consisting of 11 observables which currently set the border of the excluded part of the parameter space has been compiled in Table I. These observables have a potential to uncover the gauge LQ signal even with a small improvement of the experimental sensitivity.
For the decays P 0 → l + l ′− not listed in the catalogue, we have found the future experimental bounds needed in order to further probe the considered model. Furthermore, we have explored a class of SU (4) C models with extra heavy vector-like leptons and searched for additional possible future first signals of the gauge LQ. We have also found the smallest of these models capable of accommodating the neutral current anomalies in B decays and identified the key future measurement which can exclude such a setup. ACKNOWLEDGEMENT We acknowledge the support from the Grant agency of the Czech Republic, project no. 20-17490S, from the Grant Agency of Charles University (GAUK) project no. 12481/2019, and from the Charles University Research Center UNCE/SCI/013. We would like to express our gratitude to Michal Malinský for his valuable advice and all the support. 3. Inspecting non-diagonal parts of the anticipated approximate LF symmetry consists especially in testing the lepton flavour universality (LFU) which can be associated with the group of permutation matrices (S 3 ) LFU ⊂ U (3) LF .
Since neither U (1) 2 LF nor (S 3 ) LFU is a subgroup of the other, LFV does not necessarily imply LFU violation (LFUV) nor vice versa.
Let us trace the fate of these would-be symmetries in leptoquark interactions. For clarity of expression, consider only a single term, sayd R / U 1 V RêR ; the generalization to full-fledged interaction such as those in Eq. (2) is straightforward.
1. Apparently, the LQ interaction with the leptons and quarks conserves the lepton number L regardless of the form of the interaction matrix V R , provided the U 1 leptoquark carries L = −1.
2. If two columns of the interaction matrix V R are zero, then the LQ can be ascribed the corresponding flavour number (L e , L µ or L τ ) and there is no LFV. In the case V R has a single zero column, only a one-dimensional subgroup of U (1) 2 LF is a symmetry of the interaction (only the non-interacting flavour remains preserved). If all its columns are non-empty, U (1) 2 LF is completely explicitly broken. 3. On the other hand, respecting the (S 3 ) LFU symmetry requires that all three columns of V R are equal. Thus, the leptoquark brings new sources of LFUV whenever (at least) two columns of V R differ.
These observations hold generally, for any kind of LQ and its interaction matrix. In principle, the form of the interaction matrices may be such that either U (1) 2 LF or (S 3 ) LFU is an exact symmetry of the LQ interactions.
However, in the particular case of the gauge LQ in quark-lepton unification, the interaction matrix V R is a subject of the unitarity conditions: the column normalization rule implies that none of the columns can be empty, the U (1) 2 LF symmetry is completely broken and the LQ inevitably mediates LFV processes. Complementarily, the column orthogonality condition implies violation of (S 3 ) LFU .
In fact, no nontrivial subgroup of SU (3) LF can be a symmetry ofd R / U 1 V RêR for any invertable (e.g. unitary) V R : assuming X ∈ U (3) LF acts asê R → Xê R and U 1 → e iφ(X) U 1 , the considered interaction remains intact if and only if e iφ(X) V R X = V R , i.e., if X is a mere phase.
is considered. The details of the sequential breaking of the G 421 symmetry including this step are summarized in Table VI. In the first step of symmetry breaking, the SU (4) C factor is spontaneously broken at some high scale way above the electroweak one, which (unlike for GUTs) can be chosen arbitrarily since our framework unifies the fermions but not the gauge interactions. The smallest possible first step of the SU (4) C breaking is The Abelian factor in Eq. (B2) is generated by Quite commonly, its multiple is being used instead, which is compensated by redefinition of the gauge coupling: g BL = 3/8 g 4 . The name of the [B −L] generator is motivated by its action on the unified fermionic representations in Eq. (1). However, one must keep in mind that action of this symmetry generator in Eq. (B4) does not necessarily coincide with the difference between the baryon number B and lepton number L for other fields in the model. For example, in the Minimal Quark-Lepton Symmetry Model [3], both B and L are perturbatively conserved to all orders while T 15 C is spontaneously broken. Next, the extended models studied in Section V contain leptonic fields 1 ℓ and 1 e which transform trivially under SU (4) C and hence also under U (1) [B−L] . To emphasise the distinction between B − L and the gauge symmetry generator (B4), we shall keep the square brackets around the latter in order to indicate that [B −L] is an indivisible symbol.
The 19 gauge fields of the model can be cast as follows: In Eq. (B5a), the (3+1) × (3+1) block notation has been used. Together with the gluons G and charged intermediate vector bosons W ± , one can easily identify the vector leptoquark U 1 . Furthermore, the three electrically neutral fields A 15 , B ′ , and W 3 mix into the photon, the Z boson, and to Z ′ . The symmetry breaking (B2) gives mass only to the gauge leptoquark; the Z ′ boson acquires mass no sooner than during the second step, Thus, while the precise ratio of m U1 /m Z ′ depends on the scalar sector of the model, Z ′ can never be much heavier than U 1 . The rotation of the electrically neutral gauge fields to the mass basis can be written as where tan θ ′ = g R /(2g BL ) at the relevant scale, and θ W is the weak mixing angle (see Table VI). The angle θ m is very small when the symmetry breaking (B6) occurs way above the electroweak energy scale [3]. Hence, in the limit m Z ′ /m Z → ∞, the Z ′ boson is given by and the Z ′ coupling can be obtained by rewriting the relevant terms in the covariant derivative using relations from Table VI, which is an analogue to the SM case In the models of QLU, where all the fermions arise from SU (4) C quadruplets, the Z ′ interactions with both quarks and leptons are flavour-diagonal and universal, i.e., they respect the entire U kinematic region. NP contributions to these processes are constrained by the high-p T dilepton spectra measurements by Atlas and CMS, leading to limits around m Z ′ > 5 TeV (depending on the Z ′ coupling assumed) [86,87]. As noted in Ref. [60], these limits also indirectly constrain the mass of the gauge LQ. This bound is important in models accommodating the anomalous value of R D which require m U1 ∼ 2 TeV.
Ref. [60] further states that "the couplings of the Z ′ to SM fermions are necessarily flavour universal " and "proportional to the identity matrix in flavour space" even in the models with extra fermions because the relevant charged lepton mixing "necessarily involve states with the same B − L charge". This is, however, a misconception arising from not-distinguishing between the gauge symmetry generator [B −L] and the difference of the accidental global symmetries B − L. All the fermionic fields 4 ℓ L , 4 e R , 1 ℓ L,R , 1 e L,R are fully justified to be called leptons and carry the lepton number L, which is conserved by the gauge interactions. On the other hand, only the fields 4 ℓ L and 4 e L , which stem from SU (4) C quadruplets, are also charged with respect to [B −L], the diagonal generator of the SU (4) C group.
As a consequence of this, rotating the the left-handed [B −L] − 2Y sin 2 θ ′ lepton currents into the mass basis [see Eq. (21)] yields and similarly for the right-handed currents: Finally, using the implicit definition of V L,R in Eq. (21) and the block notation of Eq. (22), one arrives to the following formula for the Z ′ couplings with the SM fermions: where s ′2 ≡ sin 2 θ ′ amounts to 0.08 at the 2 TeV scale or to s ′2 ≃ 0.12 in the 200 TeV ballpark (assuming SMlike gauge coupling running up to m Z ′ ). Thus, the Z ′ interactions with leptons in the extended SU (4) C models are not necessarily flavour-universal and, in general, the diagonal couplings could actually be strongly suppressed. As a consequence, the limits on m Z ′ from the highenergy dilepton spectra may be considerably weakened for certain patterns of V 0 L,R . The simplified reasoning of Ref. [60] mentioned above has been used as a no-go argument for abandoning the models with the G 421 gauge group and focusing on G 4321 -based models instead when attempting to accommodate R D ( * ) . In this respect, we note that achieving the form of V 0 L,R from Eq. (4) in the framework of extended G 421 models would imply that the Z ′ couplings to the e and µ leptons are suppressed. Since the Z ′ * → τ + τ − channel is experimentally less constrained [88], a valid no-go argument needs to be more subtle. Nevertheless, the scenarios with the SU (4) Cbreaking scale as low as 2 TeV require full model specification since the effects of the new scalar and fermionic degrees of freedom would be important. This is far beyond the scope of this paper.
In any case, this study is focusing on the extended SU (4) C models with k L + k R ≤ 2. Such frameworks can not accommodate the R D ( * ) anomalies even if the Z ′ is completely ignored due to the residual constraints on the leptoquark interaction matrices V 0 L,R from the unitarity of V L,R .
During scanning of the parameter space of these models, we have not encountered a parameter point allowing for m U1 smaller than 18 TeV. Since the models allow for a similarly heavy Z ′ , the constraints from this field are not severe: unlike the gauge LQ, Z ′ does not contribute to the scalar-type 2-quark-2-lepton operators O ℓedq but only to the Wilson coefficients multiplying the vectortype ones (O ed , O ℓq ) and further to those of flavourconserving 4-lepton or 4-quark operators, all of which are experimentally less restricted.
In this analysis, the Z ′ contributions to the Wilson coefficients are not calculated. Including them could be a part of a future study focusing on the extended SU (4) C models.
Appendix C: Optimizing the scanning procedure The experimental data collected over the last decades provided rather stringent constraints on mass of the considered leptoquark. Some of the most restraining processes are the decays of K 0 L → l + l ′− and the µ → e conversion on gold nuclei, see Table I. Here we identify areas in the parameter space in which these decays are suppressed, and thus allow for lighter leptoquark.
The scanning procedure mentioned in Section III C is optimized when we restrict the parameter space to a subspace in which Schematically, the V L and V R matrices read where u, d, s are quarks, e, µ, τ are leptons and each element represents the strength of interaction of these two fermions with the leptoquark. The block matrices V I , V II , and V III are present only in the extended models studied in Section V. As follows from Eqs. (9) -(12), the BR's of the leptonic K 0 L decays are proportional to All these β's vanish if and only if .
When we think of V Lql as fixed numbers and V Rql as the unknowns, the necessary condition for a nontrivial solution to exist shrinks to where |V | stands for determinant. On the other hand, we can treat V Rql as fixed numbers and V Lql as variables, which leads to analogous result for the V R matrix. Hence, the determinants of the top left 2 × 2 submatrices of both V L and V R has to be equal to zero, regardless of the dimensionality of these matrices. Now we use a simplified rule of Laplace expansion in multiple rows (as derived by Laplace in 1772, more on this e.g. [89]) where |M | IJ is the IJ-minor, i.e., the determinant of the submatrix obtained from M by deleting rows and columns from sets I, J ⊂ D = {1, . . . , dim(M )} . The set I ′ (J ′ ) is the complement of I (J) in D, so that every row and every column index appears exactly once in Eq. (C6). If M is a unitary matrix, its determinant is just a complex phase, and we can further simplify Eq. (C6) to In other words, the determinant of any submatrix of a unitary matrix is equal in magnitude to the determinant of the complementary submatrix. Applying this observation to the V L and V R matrices of dimension 3, Eq. (C5) leads to with an implication BR V (B 0 d,s → τ τ ) = 0 (which has been derived also in Ref. [4] by directly solving Eqs. (C3) in a specific parametrization). Note also that the "solutions" to the B anomalies which leave out the unitarity constrains such as Eq. (4), usually require that the V L,R bτ elements are the largest ones in order to address also R D ( * ) . Therefore, the anomalous value of R D ( * ) can not be accounted to the Pati-Salam-type leptoquark.
Fulfilling the rather simple condition (C5) can be tough for V L , V R of higher dimensions. To this end, we introduce the composite parametrization of U (n) matrices [90,91], which turns out to be particularly convenient in this respect. Its n 2 parameters λ ij consist of 1 2 n(n − 1) angles (i < j) and 1 2 n(n + 1) phases (i ≥ j). where c ij = cos λ ij and s ij = sin λ ij . For higher dimensions we refer to the original literature [90,91] or to the publicly available implementation in Wolfram Mathematica [92]. The slice of the parameter space used for Fig. 3 is obtained by fixing these parameters according to Table VII.
The beautiful advantage of the composite parametrization is that the subspace obeying Eq. (C1) can be obtained by fixing the same parameters for any dimension ≥ 3 of V L and V R . The necessary condition (C5) for existence of the solution is fulfilled by With this in hand, it can be shown that setting solves Eqs. (C4) entirely. Naïvely, other solutions can be found but they fall outside the proper domain of the λ's. An equivalent solution to (C3) was found in Ref. [4] for dim V L = dim V R = 3 within a different parametrization.
2. Avoiding CR(µ → e, Au) Leaving the K 0 L decays for a while, we now focus on another very important constraint stemming from the limits λ11L λ12L λ13L λ21L λ22L λ23L λ31L λ32L λ33L TABLE VII. Angles and phases defining the slice of the parameter space used in Fig. 3. on µ → e conversion on gold nuclei, CR(µ → e, Au) < 7 × 10 −13 [43]. In the same manner, we enforce CR(µ → e, Au) = 0. (C12) A leptoquark with Q = + 2 /3 mediates this process at the tree level by an interaction with the d quarks and the sea s quarks in the nucleons. The calculation in flavio is based on Ref. [93]. The scalar-type effective vertices, (d R d L )(e L µ R ) and (d R d L )(µ L e R ), are predicted to engage in this process even more efficiently than the vector-type ones. Thus, to avoid these constraints when searching for limits from other interesting processes, the following condition must be approximately fulfilled: It can be shown that any V L,R pair obeying Eq. (C13) together with the set of Eqs. (C3) must necessarily have some of the elements from the upper left 2 × 2 submatrix equal to zero. The possible patterns for V L,R are where • denotes an unfixed value. The last two cases are available only when V L or V R has dimension n ≥ 4, respectively.
Finding the unitary parametrization fulfilling both Eqs. (C3) and (C13) is straightforward though somewhat tedious as the solution has to be found for each dimension of V L and V R separately.
Notable but order-of-magnitude smaller contributions to the coherent µ → e conversion still arise from vectortype operators (triggered by V Lde V L * dµ and V Rde V R * dµ ) and well as the muon conversion on the sea s-quarks in the nucleons (such amplitudes are proportional to V Lse V R * sµ or V Lse V R * sµ ).