Addressing the RD ( ) anomalies with an S 1 leptoquark from SO ( 10 ) grand unification

Ufuk Aydemir , Tanumoy Mandal , and Subhadip Mitra 5,‡ Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, P. R. China School of Physics, Huazhong University of Science and Technology, Wuhan, Hubei 430074, P. R. China Department of Physics and Astronomy, Uppsala University, Box 516, SE-751 20 Uppsala, Sweden Department of Physics and Astrophysics, University of Delhi, Delhi 110 007, India Center for Computational Natural Sciences and Bioinformatics, International Institute of Information Technology, Hyderabad 500 032, India


I. INTRODUCTION
During the past few years, several disagreements between experiments and the Standard Model (SM) predictions in the rare B decays have been reported by the BABAR [1,2], LHCb [3][4][5][6][7] and Belle [8][9][10][11] collaborations. So far, these anomalies have been quite persistent. The most significant ones have been observed in the R D ðÃÞ and R K ðÃÞ observables, defined as,

R D ðÃÞ ¼
BRðB → D ðÃÞ τνÞ BRðB → D ðÃÞ lνÞ and Here, l ¼ fe or μg and BR stands for branching ratio. The experimental values of R D and R D Ã are in excess of their SM predictions [12][13][14][15] by 1.4σ and 2.5σ, respectively, (combined excess of 3.08σ in R D ðÃÞ ) based on the world averages as of spring 2019, according to the Heavy Flavor Averaging Group [16], whereas, the observed R K and R K Ã are both suppressed compared to their SM predictions [17,18] by ∼2.6σ. One of the possible explanations for the B-decay anomalies is the existence of scalar leptoquarks whose masses are in the few-TeV range ). 1 Leptoquarks, which posses both lepton and quark couplings, often exist in the grand unified theories (GUTs) or in the Pati-Salamtype models. Considering that the LHC searches, except for these anomalies, have so far returned empty handed, if the leptoquarks indeed turn out to be behind these anomalies, it is likely that there will only be a small number of particles to be discovered. However, their existence in small numbers at the TeV scale would be curious in terms of its implications regarding physics beyond the Standard Model. Scalars in these models mostly come in large multiplets and it would be peculiar that only one or a few of the components become light at the TeV scale while others remain heavy. Mass splitting is actually a well-speculated subject in the literature in the context of the infamous doublet-triplet splitting problem in supersymmetric theories and GUTs. From this point of view, even the SM Higgs in a GUT framework is troublesome if it turns out to be the only scalar at the electroweak (EW) scale.
In this paper, we consider a single scalar leptoquark, S 1 ð3; 1; −1=3Þ, at the TeV scale, in the SO(10) GUT framework [62][63][64][65][66][67][68][69][70][71][72][73][74][75]. This particular leptoquark was discussed in the literature to be responsible for one (or possibly both) of the R D ðÃÞ and R K ðÃÞ anomalies [21,22,25,32]. (Later, however, it was shown in Ref. [45] that a single S 1 leptoquark can only alleviate the R K ðÃÞ discrepancy, but cannot fully resolve it.) Furthermore, it is contained in a relatively small multiplet that resides in the fundamental representation of the SO(10) group, a real 10, together with a scalar doublet with the quantum numbers that allow it to be identified as the SM Higgs. Therefore, S 1 being the only light scalar entity other than the SM Higgs doublet is justified in this scenario. A future discovery of such a leptoquark at the TeV scale could be interpreted as evidence in favor of an SO(10) GUT.
It has been argued in the literature that a real 10 H in the minimal SO(10) setup (even together with a real 120 H or a complex 120 H ) is not favored in terms of a realistic Yukawa sector [68]. On the other hand, it has been recently discussed in Ref. [75] that a Yukawa sector consisting of a real 10 H , a real 120 H and a complex 126 H , can establish a realistic Yukawa sector due to the contributions from the scalars whose quantum numbers are the same as the SM Higgs doublet. Thus, this is the scalar content we assign in our model for the Yukawa sector.
The inclusion of S 1 in the particle content of the model at the TeV scale does not improve the status of the SM in terms of gauge coupling unification, which cannot be realized by the particle content in question. Fortunately, in the SO(10) framework, there are other ways to unify the gauge coupling constants, in contrast to models based on the SU(5) group which also contains such a leptoquark within the same multiplet as the SM Higgs. As we illustrate in this paper, inserting a single intermediate phase where the active gauge group is the Pati-Salam group, SUð4Þ C ⊗ SUð2Þ L ⊗ SUð2Þ R , which appears to be the favored route of symmetry breaking by various phenomenological bounds [71], establishes coupling unification as desired. In our model, the Pati-Salam group is broken into the SM gauge group at an intermediate energy scale M C , while SO (10) is broken into the Pati-Salam group at the unification scale M U . We consider two versions of this scenario, depending on whether the left-right symmetry, so-called D parity, is broken together with SO(10) at M U , or it is broken at a later stage, at M C , where the Pati-Salam symmetry is broken into the SM gauge symmetry.
Light color triplets, similar to the one we consider in this paper, are often dismissed for the sake of proton stability since these particles in general have the right quantum numbers for them to couple potentially dangerous operators that mediate proton decay. On the other hand, the proton stability could possibly be ensured through various symmetry mechanisms such as the utilization of Peccei-Quinn (PQ) symmetry [68,76], other U(1) symmetries such as the one discussed in Ref. [77], or a discrete symmetry similar to the one considered in Ref. [22]. Operators leading to proton decay could also be suppressed by a specific mechanism such as the one discussed in Ref. [78] or they could be completely forbidden by geometrical reasons [38]. In this paper, we adopt a discrete symmetry as suggested in Ref. [22], assumed to operate below the intermediate symmetry-breaking scale even though it is not manifest at higher energies.
Motivated by the possible existence of a single TeV-scale S 1 in the SO(10) GUT framework, we move on to investigate the phenomenological implications of our model. In Ref. [22], it was shown that a TeV-scale S 1 leptoquark can explain the R D ðÃÞ anomalies while simultaneously inducing the desired suppression in R K ðÃÞ through box diagrams. Since the most significant anomalies are seen in the R D ðÃÞ observables, in this paper, we concentrate mostly on scenarios that can accommodate these observables. Generally, a TeV-scale S 1 requires one Yukawa coupling to be large to accommodate the R D ðÃÞ anomalies [32,45]. This, however, could create a problem for the b → sνν transition rate measured in the R νν K observable. In the SM, this decay proceeds through a loop whereas S 1 can contribute at the tree level in this transition. Therefore, the measurement of R νν K is very important to restrict the parameter space of S 1 . 2 Some specific Yukawa couplings of S 1 are also severely constrained from the Z → ττ decay [47] and the LHC ττ resonance search data [48]. Therefore, it is evident that in order to find the R D ðÃÞ -favored parameter space while successfully accommodating other relevant constraints, one has to introduce new d.o.f. in terms of new couplings and/or new particles. Here, we consider a specific Yukawa texture with three free couplings to show that a TeV-scale S 1 , consistent with relevant measurements, can still explain the R D ðÃÞ anomalies.
The rest of the paper is organized as follows. In Sec. II, we introduce our model. In Sec. III, we display the unification of the couplings for two versions of our model. In Sec. IV, we present the related LHC phenomenology with a single extra leptoquark S 1 . We display the exclusion limits from the LHC data and discuss related future prospects. We also study the renormalization group (RG) running of the Yukawa couplings. Finally in Sec. V, we end our paper with a discussion and conclusions.

II. THE SO(10) MODEL
In our SO(10) model, we entertain the idea that the SM Higgs doublet is not the only scalar multiplet at the TeV 2 One can avoid this conflict by introducing some additional degrees of freedom (d.o.f.), as shown in Ref. [30]. There, the authors introduced an S 3 leptoquark in addition to the S 1 to concomitantly explain R D ðÃÞ and R K ðÃÞ while being consistent with R νν K .
scale, but it is accompanied by a leptoquark S 1 ¼ ð3; 1; −1=3Þ, both of which reside in a real ten-dimensional representation, 10, of SO(10) group. The peculiar mass splitting among the components of this multiplet does not occur, leading to a naturally light scalar leptoquark at the TeV scale. We start with a real 10 H of SO(10) whose Pati-Salam and SM decompositions are given as 10 ¼ ð1; 2; 2Þ 422 ⊕ ð6; 1; 1Þ 422 where subscripts denote the corresponding gauge group and we set The scalar content we assign for the Yukawa sector consists of a real 10 H , a real 120 H , and a complex 126 H , which establishes a realistic Yukawa sector through mixing between the scalars whose quantum numbers are the same as the SM Higgs doublet, as shown in Ref. [75]. Note that where 16 is the spinor representation in which each family of fermions, including the right-handed neutrino, resides in. The subscripts s and a denote the symmetric and antisymmetric components. The Pati-Salam decompositions of 120 and 126 are given as The Yukawa terms are then given as where Y 10 and Y 126 are complex Yukawa matrices, symmetric in the generation space, and Y 120 is a complex antisymmetric one. The SM doublet contained in ϕð1; 2; 2Þ 422 of 10 is mixed with other doublets accommodated in ϕð1; 2; 2Þ 422 and Σð15; 2; 2Þ 422 of the real multiplet 120 and Σð15; 2; 2Þ 422 of 126, yielding a Yukawa sector consistent with the observed fermion masses [75]. The fermion mass matrices for the up quark, down quark, charged leptons, Dirac neutrinos and Majorana neutrinos are given as [75] 120 ÞY 120 ; 120 ÞY 120 ; 120 ÞY 120 ; where v u i ðv d i Þ are the presumed vacuum expectation values of the neutral components of the corresponding bidoublets, contributing to the masses of the up quarks and the Dirac neutrinos (the down quarks and the charged leptons), and where , and v L ∼ v R v 2 10 =M 2 GUT is induced by the potential term 10 2 H 120 2 H and is responsible for the left-handed Majorana neutrino masses (type-II seesaw) [75].
For the symmetry-breaking pattern, we consider a scenario in which the intermediate phase has the gauge symmetry of the Pati-Salam group SUð4Þ C ⊗ SUð2Þ L ⊗ SUð2Þ R . Note that this specific route of symmetry breaking appears to be favored by various phenomenological bounds [71]. We consider two versions of this scenario depending on whether or not the Pati-Salam gauge symmetry is accompanied by the D-parity invariance, a Z 2 symmetry that maintains the complete equivalence of the left and right sectors [65,66,79], after the SO(10) breaking. The symmetry-breaking sequence is schematically given as where, we use the notation, In the second scenario, which we call model A 2 , the first stage of the symmetry breaking is realized through the Pati-Salam singlet contained in 54, which acquires a VEV. This singlet is even under D parity, and therefore, D parity is not broken at this stage with the SO(10) symmetry, and the resulting symmetry group valid down to M C is G 422 . The rest of the symmetry breaking continues in the same way as in model A 1 . Consequently in model A 2 , we include one more Pati-Salam multiplet at M C , Δ L ð10; 3; 1Þ, in order to maintain a complete left-right symmetry down to M C . The scalar content and, for later use, the corresponding RG coefficients in each energy interval are given in Table I.
Finally, the relevant Lagrangian for phenomenological analysis at low energy is given by where Q and L are the SM quark and lepton doublets (for each family), Λ L=R are coupling matrices in flavor space and ψ c ¼ Cψ T are charge-conjugate spinors. Notice the absence of dangerous diquark couplings of S 1 that would lead to proton decay. One way to forbid these couplings is to impose a Z 2 symmetry [22] that emerges below the Pati-Salam breaking scale under which quarks and leptons transform with opposite parities whereas the leptoquark is assigned odd parity, i.e., ðq; l; S 1 Þ → ðAEq; ∓l; −S 1 Þ. Note also that the inclusion of S 1 can affect the stability of the electroweak vacuum via loop effects. The relevant discussion can be found in Ref. [80]. Evidently, the Lagrangian given in Eq. (9) together with the SM Lagrangian should be understood in the effective field theory context. The new and SM Yukawa couplings in the TeV-scale Lagrangian are induced from the original SO (10) Yukawa couplings each of which is generated by a linear combination of unification-scale operators and gets modified due to the mixing effects induced by the scalar fields that have Yukawa couplings to three chiral families of 16 F . It is indeed this rich structure that enables the realization of a fermion mass spectrum consistent with the expected fermion masses of the SM model at the unification scale, as shown in Ref. [75]. As we will discuss later in Sec. IV G, the modification to the SM RG running of the Yukawa couplings due to the inclusion of Λ L=R does not register strong changes in the fermionic mass spectrum and hence the main message of Ref. [75] is valid in our case, as well.

III. GAUGE COUPLING UNIFICATION
In this section, after we lay out the preliminaries for oneloop RG running and show that the new particle content at the TeV scale does not lead to the unification of the SM gauge couplings directly, we illustrate gauge coupling unification with a single intermediate step of symmetry breaking. Once the particle content at low energies is determined, there may be numerous ways to unify the gauge couplings, depending on the selection of the scalar content in SO(10) representations. In the literature, the canonical way to make this selection is through adopting a minimalistic approach, allowed by the observational constraints. In the following, we pursue the same strategy while taking into account the analysis made in Ref. [75] for a realistic Yukawa sector.
A. One-loop RG running where the RG coefficients a i are given by [81,82] and the full gauge group is given as The summation in Eq. (11) is over irreducible chiral representations of fermions (R f ) and irreducible representations of scalars (R s ) in the second and third terms, respectively. The coefficient η is either 1 or 1=2, depending on whether the corresponding representation is complex or (pseudo)real, respectively. d j ðRÞ is the dimension of the representation R under the group G j≠i . C 2 ðG i Þ is the quadratic Casimir for the adjoint representation of the group G i , and T i is the Dynkin index of each representation (see Table II). For the U(1) group, C 2 ðGÞ ¼ 0 and where Y is the Uð1Þ Y charge. The addition of S 1 to the particle content of the SM does not help in unifying the gauge couplings as displayed in Fig. 1, where the RG running is performed with the modified RG coefficients given in Table I

B. Unification with a single intermediate scale
We start by labeling the energy intervals in between symmetry-breaking scales ½M Z ; M C and ½M C ; M U with Roman numerals as I∶ ½M Z ; M C ; G 213 ðSMÞ; The boundary/matching conditions we impose on the couplings at the symmetry-breaking scales are   (3), and SU(4). Our normalization convention in this paper follows the one adopted in Ref. [82]. Notice that there are two inequivalent 15-dimensional irreducible representations for SU (3). We use the central values of the low-energy data as the boundary conditions in the RG running (in the MS scheme) [83,84]: The coupling constants are all required to remain in the perturbative regime during the evolution from M Z to M U . The RG coefficients, a i , differ depending on the particle content in each energy interval, changing every time symmetry breaking occurs. Together with the matching and boundary conditions, one-loop RG running leads to the following conditions on the symmetry-breaking scales M U and M C : where the notation on a i is self-evident. The unified gauge coupling α U at the scale M U is then obtained from Thus, once the RG coefficients in each interval are specified, the scales M U and M C , and the value of α U are uniquely determined. The results are given in Table III, and unification of the couplings is displayed in Fig. 2 for each model. As mentioned previously, we assume in this paper that the proton-decay-mediating couplings of S 1 are suppressed. On the other hand, we do not make any assumptions regarding the other potentially dangerous operators which could lead to proton decay. Thus, it is necessary to inspect whether the predictions of our models displayed in Table III are compatible with the current bounds coming from the proton decay searches or not. The most recent and stringent bound on the lifetime of the proton comes from the mode p → e þ π 0 , and is τ p > 1.6 × 10 34 years [85]. As for the proton decay modes that are mediated by the super-heavy gauge bosons, which reside in the adjoint representation of SO (10) , we obtain M U ≳ 10 15.9 GeV, which is consistent with predictions of both model A 1 and model A 2 , within an order of magnitude of the latter. Additionally, since M C is the scale at which the Pati-Salam symmetry breaks into the SM, it determines the expected mass values for the proton-decaymediating color triplets. From a naive analysis [71], it can be shown that the current bounds on the proton lifetime require M C ≳ 10 11 GeV, again consistent with the predictions of both model A 1 and model A 2 , within an order of magnitude of the former. Note that these bounds should be taken as order-of-magnitude estimates since, while obtaining them, we approximate the anticipated masses of the super-heavy gauge bosons and the color triplets as M X ≈ M U and M T ≈ M C , while it would not be unreasonable to expect that these mass values could differ from the corresponding energy scales within an order of magnitude.

IV. LOW-ENERGY PHENOMENOLOGY
The existence of a TeV-scale charge −1=3 scalar leptoquark in a GUT framework is quite interesting from a phenomenological perspective mainly for two reasons. First, its existence is testable at the LHC. The directdetection searches for scalar leptoquarks have been putting exclusion bounds on S 1 with different decay hypotheses [87,88]. Second, as mentioned earlier, such a leptoquark can offer an explanation of some persistent flavor anomalies observed in several experiments. For example, if we consider the anomalies observed in the B-meson semileptonic decays via charged currents (collectively these show the most significant departure from the SM expectations), S 1 can provide an explanation if it couples with τ and neutrino(s) and b and c quarks. The direct LHC bounds on such a leptoquark are not very severe but as it has been pointed out in Ref. [48], the present LHC data in the pp → ττ=τν channels have actually put constraints on the S 1 parameter space relevant for explaining the observed R D ðÃÞ anomalies. Here, using flavor data and LHC constraints, we obtain the allowed parameter space in our model. We also point out some possible new search channels at the LHC. On the flavor side, our primary focus is on the chargedcurrent anomalies observed in the semileptonic B decays in the R D ðÃÞ observables. Hence, as in Ref. [48], we focus on the interaction terms of S 1 that could play a role to address the R D ðÃÞ anomalies for simplicity.

A. The S 1 model
The single TeV-scale S 1 leptoquark that originates from the GUT model discussed in Sec. II transforms under the SM gauge group as ð3; 1; −1=3Þ. The low-energy interactions of S 1 with the SM fields are shown in a compact manner in Eq. (9). Below, we display the relevant interaction terms required for our phenomenological analysis, where Q i and L i denote the ith-generation quark and lepton doublets, respectively and λ H ij represents the coupling of S 1 with a charge-conjugate quark of the ith generation and a lepton of the jth generation with chirality H. Without any loss of generality, we assume all λ's are real in our collider analysis since the LHC data that we consider are insensitive to their complex nature. Also, we only consider mixing among quarks [Cabibbo-Kobayashi-Maskawa (CKM) mixing] and ignore neutrino mixing (Pontecorvo-Maki-Nakagawa-Sakata mixing) completely as all neutrino flavors contribute to the missing energy and hence are not distinguishable at the LHC. The couplings of S 1 to the first-generation SM fermions are heavily constrained [32]. Hence, we assume λ 1i ; λ i1 ¼ 0 in our analysis. 3 The parton-level Feynman diagrams for the b → cτν decay (responsible for the B → D ðÃÞ τν decay) are shown in Fig. 3. In order to have a nonzero contribution in the R D ðÃÞ observables from S 1 , we need the bνS 1 and cτS 1 couplings to be nonzero simultaneously. Minimally, one can start with just a single free coupling: either λ L 23 or λ L 33 . The coupling λ L 23 (λ L 33 ) directly generates the cτS 1 (bνS 1 ) interaction and the other one, i.e., the bνS 1 (cτS 1 ) can be generated through the CKM mixing among quarks. These two minimal scenarios were discussed in detail in Ref. [48]. For these two cases, the Lagrangian in Eq. (17) can be written explicitly as, In Ref. [48], it was shown that for C 1 , the R D ðÃÞ -favored parameter space is already ruled out by the latest LHC data. On the other hand, C 2 is not seriously constrained by the LHC data since this scenario is insensitive to the coupling λ L 33 . Only the pair-production searches, which are largely The parton-level process in the presence of S 1 that would contribute to B → Kνν is shown in diagram (c). 3 However, these couplings can be generated through the CKM mixing. We refer the interested readers to Ref. [32] for various important flavor constraints in this regard. insensitive to λ L 33 , in the ttττ and bbνν modes exclude M S 1 up to 900 GeV [87] and 1100 GeV [88], respectively for a 100% BR in each decay mode. However, Ref. [47] showed that the R D ðÃÞ -favored parameter space in C 2 is also ruled out by the electroweak precision data on the Z → ττ decay.
The above two minimal cases, C 1 and C 2 , are the two extremes. One can, however, consider a next-to-minimal situation, where both λ L 23 and λ L 33 are nonzero to explain R D ðÃÞ anomalies being within the LHC bounds [48]. However, B → K ðÃÞ νν decay results severely constrain such a scenario due to the tree-level leptoquark contribution [see Fig. 3(c)]. References [32,45] indicated that a large λ R 23 might help explain various flavor anomalies simultaneously while being consistent with other relevant experimental results.
In this paper, we allow λ L 23 , λ L 33 and λ R 23 to be nonzero and perform a parameter scan for a single S 1 solution of the R D ðÃÞ anomalies. We locate the R D ðÃÞ -favored parameter space that satisfies the limits from B → K ðÃÞ νν and Z → ττ decays and is still allowed by the latest LHC data. A S 1 can also provide new final states at the LHC like ττ þ jets and τ þ = E T þ jets in which leptoquarks have not been searched for before.
In the SM, the semitauonic B decay is mediated by the left-handed charged currents and the corresponding four-Fermi interactions are given by the following effective Lagrangian: In the presence of new physics, there are a total of five four-Fermi operators that appear in the effective Lagrangian for the B → D ðÃÞ τν decay [89], where the C X 's are the Wilson coefficients associated with the effective operators: (1) Vector operators: (2) Scalar operators: (3) Tensor operator: The operator O V L is SM-like and the other four operators introduce new Lorentz structures into the Lagrangian. Note that the operator O T R is identically zero, i.e., The S 1 leptoquark that we consider can generate only O V L ;S L ;T L . Hence, the coefficients of the other two operators, namely, C V R and C S R remain zero in our model. In terms of the S 1 parameters the Wilson coefficients can be expressed as, These relations are obtained at the mass scale M S 1 . However, running of the strong coupling constant down to m b ∼ 4.2 GeV changes these coefficients substantially except for C V L which is protected by the QCD Ward identity. As a result, the ratio C S L =C T L becomes, The modification factor ρ can be obtained from Ref. [32], and we display it in Fig. 4. In terms of the nonzero Wilson coefficients we can express the ratios r D ðÃÞ ¼ R D ðÃÞ =R SM With Eq. (24) one can simplify the above equations as, where ρ ¼ ρðm b ; M S 1 Þ. There are two other observables related to the R D Ã -the longitudinal D Ã polarization F L ðD Ã Þ and the longitudinal τ polarization asymmetry P τ ðD Ã Þ-which have recently been measured by the Belle Collaboration [10,11,90]. In terms of the nonzero Wilson coefficients in our model, F L ðD Ã Þ and P τ ðD Ã Þ are expressed as [49], These equations can further be simplified as, These two observables have the power to discriminate between new physics models with different Lorentz structures (see e.g., Ref. [91]). In Table IV, we list the bounds on the R D ðÃÞ related observables that we include in our parameter scan.

C. Constraint from R νν
The SM flavor-changing neutral-current b → sνν transition proceeds through a loop and is suppressed by the Glashow-Iliopoulos-Maiani mechanism, whereas in our model, S 1 can mediate this transition at the tree level [see Fig. 3(c)]. Therefore, this neutral-current decay can heavily constrain the parameter space of our model. We define the following ratio: The current experimental 90% confidence limit (C.L.) upper limits on the above quantities are R νν K < 3.9 and R νν K Ã < 2.7 [93]. In terms of our model parameters, R νν K ðÃÞ is given by the following expression: where a ¼ ffiffi ffi 2 p π 2 =ðe 2 G F jC SM L jÞ with C SM L ≈ −6.38 [32]. We use this constraint in our analysis and find that it significantly restricts our parameter space. Note that this constraint applies on λ L 23 and λ L 33 but not on λ R 23 . In Fig. 5, we show the regions in the λ L 23 − λ L 33 plane with R νν K Ã < 2.7 for two different values of M S 1 .

D. Constraint from Z → ττ decay
Another important constraint comes from the Zττ coupling measurements. The Z → ττ decay is affected by the S 1 loops as shown in Ref. [47]. The contribution of S 1 to the Zττ coupling shift (Δκ Zττ ) comes from a loop with an up-type quark (q) and an S 1 . The shift scales as the square of the S 1 tτ coupling (λ L=R q3 ) and m 2 q . Hence, the dominant contribution comes from when q is the top quark implying that the Zττ coupling measurements can restrict only λ L 33 but not λ L 23 or λ R 23 . For instance, we see from Ref. [47] that λ L 33 ≳ 1.4 can be excluded for M S 1 ∼ 1 TeV with 2σ confidence. We incorporate this bound into our parameter scan.

E. LHC phenomenology and constraints
We now make a quick survey of the relevant LHC phenomenology of a TeV-range S 1 that couples with τ, ν and s and c quarks. For this discussion we compute all the necessary cross sections using the universal FeynRules output [94] model files from Ref. [48] in MADGRAPH5 [95]. We use the NNPDF23LO [96] parton distribution functions. Wherever required, we include the next-to-leadingorder QCD K-factor of ∼1.3 for the pair production in our analysis [97].

Decay modes of S 1
For nonzero λ L 23 , λ L 33 and λ R 23 , S 1 can decay to cτ, sν, tτ and bν states. CKM mixing among quarks enables decays to uτ and dν but we neglect them in our analysis as the offdiagonal CKM elements are small. The BRs of S 1 to various decay modes vary depending on the coupling strengths. If λ R 23 ≫ λ L 23 ; λ L 33 , the dominant decay mode is S 1 → cτ, whereas for λ L 23 ≫ λ R 23 ; λ L 33 , BRðS 1 → cτÞ ≈ BRðS 1 → sνÞ ≈ 50%. On the other hand, when λ L 33 ≫ λ L 23 ; λ R 23 , the dominant decay modes are S 1 → tτ and S 1 → bν with about a 50% BR in each mode. Since partial decay widths depend linearly on M S 1 , BRs are insensitive to the mass of S 1 .

Production of S 1
At the LHC, S 1 can be produced resonantly in pairs or singly and nonresonantly through indirect production (t-channel S 1 exchange process).
Pair production: The pair production of S 1 is dominated by the strong coupling and, therefore, it is almost model independent. The mild model dependence enters in the pair production through the t-channel lepton or neutrino exchange processes. However, the amplitudes of those diagrams are proportional to λ 2 and generally suppressed for small λ values (for bigger M S 1 and large λ values, this part could be comparable to the model-independent part of the pair production). Pair production is heavily phase-space suppressed for large M S 1 and we find that its contribution is very small in our recast analysis. Pair production can be categorized into two types depending on the final states: symmetric, where both leptoquarks decay to the same modes, and asymmetric, where the two leptoquarks decay via two different modes. These two types give rise to various novel final states.
Symmetric modes: Asymmetric modes: where the curved connection over a pair of particles indicates that the pair is coming from a decay of S 1 . Searches for leptoquarks in some of the symmetric modes were already done at the LHC [87,88]. Leptoquark searches in some of the symmetric and most of the asymmetric modes are yet to be performed at the LHC. Single production: The single productions of S 1 , where S 1 is produced in association with a SM particle, are fully model dependent as they depend on the leptoquark-quark-lepton couplings. These are important production modes for large couplings and heavier masses (since single productions receive less phase-space suppression than the pair production). Depending on the final states, single productions can be categorized as follows.
Symmetric modes: Asymmetric modes: Here j stands for an untagged jet and "jets" means any number (≥1) of untagged jets. These extra jets can be either radiation or hard (genuine three-body single production processes can have sizeable cross sections; see Refs. [98][99][100] for how one can systematically compute them). As single production is model dependent, the relative strengths of these modes depend on the relative strengths of the coupling involved in the production as well as the BR of the decay mode involved. Indirect production: Indirect production is the nonresonant process where a leptoquark is exchanged in the t channel. With leptoquark couplings to τ and ν, this basically gives rise to three possible final states: ττðττ þ jetsÞ, τνð= E T þ τ þ jetsÞ and ννð= E T þ jetsÞ. The amplitudes of these processes are proportional to λ 2 . So the cross section grows as λ 4 . Hence, for an order-one λ, indirect production has a larger cross section than other production processes for large M S 1 (see Fig. 2 of Ref. [48]). However, the indirect production substantially interferes with the SM background process pp → V ðÃÞ → ll (l ¼ τ=ν). Though the interference is Oðλ 2 Þ, its contribution can be significant for a TeV-scale S 1 because of the large SM contribution (larger than both the direct production modes and the λ 4 indirect contribution, assuming λ ≳ 1). In general, the interference could be either constructive or destructive depending on the nature of the leptoquark species and its mass [101]. For S 1 , we find that the interference is destructive in nature [48]. Hence, for a TeV-scale S 1 if λ is large, this destructive interference becomes its dominant signature in the leptonic final states.

Constraints from the LHC
The mass exclusion limits from the pair-production searches for S 1 at the LHC are as follows. Assuming a 100% BR in the S → tτ mode, a recent search at the CMS detector has excluded masses below 900 GeV [87]. Similarly, for a leptoquark that decays exclusively to bν or sνð≡jν) final states, the exclusion limits are at 1100 and 980 GeV [88], respectively. However, going beyond simple mass exclusions, we make use of the analysis done in Ref. [48] for the LHC constraints. It contains the independent LHC limits on the three couplings shown in Eq. (17) as functions of M S 1 as well as a summery of the direct-detection exclusion limits.
Apart from the processes with = E T þ jets final states, all other production processes can have either ττ þ jets final states or = E T þ τ þ jets final states. Hence, the latest pp → Z 0 → ττ and pp → W 0 → τν searches at the ATLAS detector [102,103] were used to derive the constraints in Ref. [48]. There we notice that the limits on λ L 23 from the τν data are weaker that the ones obtained from the ττ data. The ττ data also constrain λ R 23 . From the earlier discussion, it is clear that the interference contribution plays the dominant role in determining these limits. However, its destructive nature means that in the signal region one would expect less events than in the SM-only predictions. Hence, the limits were obtained assuming that either λ L 23 or λ R 23 is nonzero at a time or by performing a χ 2 test of the transverse mass (m T ) distributions of the data. As, for heavy S 1 , the limits on λ L 23 and λ R 23 are dominantly determined by the interference of the indirect production, they are very similar. We can translate these limits from the ττ data on any combination of λ L 23 and λ R 23 in a simple manner assuming λ L=R . In Fig. 6 we display the limits on λ L=R 23 as a function of M S 1 . 4 The LHC data is insensitive to λ L 33 as it was shown in Ref. [48]. This can be understood from the following argument. First, the pair production is insensitive to this coupling as we have already mentioned. Second, the singleproduction process pp → tτν via an S 1 has too small a cross section (∼2 fb for λ L 33 ∼ 1 and M S 1 ¼ 1 TeV) to make any difference at the present luminosity. Finally, there is no interference contribution in the ττ and τ þ = E T channels as there is no t quark in the initial state. Hence, λ L 33 remains unbounded from these searches.

F. Parameter scan
To find the R D ðÃÞ -favored regions in the S 1 parameter space that are not in conflict with the limits on F L ðD Ã Þ, P τ ðD Ã Þ, R νν K ðÃÞ , Z → ττ decay and the bounds from the LHC, we consider two benchmark leptoquark masses: M S 1 ¼ 1 and 2 TeV. We allow all three free couplings, λ L 23 , λ R 23 and λ L 33 to vary. For every benchmark mass, we perform a random scan over the three couplings in the perturbative range − ffiffiffiffiffi ffi 4π p to ffiffiffiffiffi ffi 4π p (i.e., jλj 2 =4π ≤ 1). We do not consider complex values for the couplings. In Fig. 7, we show the outcome of our scan with different two-dimensional projections. In every plot we show two couplings and allow the third coupling to vary. In each of these plots we show the following.
(1) The flavor ∪ EW (FEW) regions: The orange dots mark the regions favored by the R D ðÃÞ observables within 95% C.L. while satisfying the available bounds on the F L ðD Ã Þ, P τ ðD Ã Þ and R νν K ðÃÞ observables (flavor bounds). In addition, these points also satisfy the bound on λ L 33 coming from the Z → ττ decay within 95% C.L. [47] (electroweak bound).
(2) The flavor ∪ EW ∪ LHC (FEWL) regions: As we take into account the limits on λ L=R 23 from the ATLAS pp → ττ data from the 13 TeV LHC [48] along with the previous constraints we obtain the regions marked by the green points. These are the points that survive all the limits considered in this paper. From the plots we see that substantial portions of parameter regions survive after all the constraints. This implies that the S 1 model can successfully explain the R D ðÃÞ anomalies. If one looks only at the R D ðÃÞ anomalies, in principle, one can just set λ L 23 and/or λ L 33 to be large. But coupling values that make C V L [see Eq. (23)] big come into conflict with the R νν K ðÃÞ bound [see Eq. (34)]. This is why we do not see any point where both λ L 23 and λ L 33 are large in the first column of Fig. 7. In addition, the LHC puts bounds on λ L=R 23 [48] whereas the Z → ττ data puts a complimentary bound on λ L 33 [47]. The restriction on λ L 33 from the Z → ττ data can be seen in Figs we see that these two couplings take opposite signs mainly because both R D Ã and F L ðD Ã Þ prefer a positive C S L [see Eqs. (28) and (31)].

G. RG running of the Yukawa couplings and perturbativity
One of the questions raised by the introduction of new Yukawa couplings is whether the new model remains perturbative up to high-enough energies. This is particularly important in the GUT framework since the RG running of the gauge couplings is performed under the assumption of perturbativity. Fortunately, with the latest LHC data, we have quite a large available parameter space in the deep perturbative region, as can be seen in Fig. 7. While this suggests that the model is safe in terms of perturbativity as long as the new Yukawa couplings are small enough at the electroweak scale, it is still informative to investigate the RG running in detail, especially the case FIG. 6. Two-sigma exclusion limits on λ L=R as a function of M S 1 as obtained in Ref. [48] from the ATLAS pp → ττ data [102]. The colored region is excluded. 4 Actually, for M S 1 between 1 and 2 TeV, the limits on λ L 23 are slightly stronger than those on λ R 23 because in the SM the Z boson couples differently to left-and right-handed τ's. However, we ignore this minor difference and take the stronger limits on λ L 23 as the limits on λ L=R 23 to remain conservative.
in which the Yukawa couplings take larger initial values. We address this issue in this part of the paper. 5 The equations for the RG running of the Yukawa couplings are given in the Appendix. Results for various benchmark cases are displayed in Fig. 8. Among the S 1 mass values we have considered in this paper, the most constrained parameter space is that of M S 1 ¼ 2 TeV. Therefore, we choose our benchmark values from the parameter space for this mass value, displayed in Fig. 7. Note that λ L 33 ðM Z Þ and λ L 23 ðM Z Þ cannot both be large and that λ L 33 ðM Z Þ and λ R 23 ðM Z Þ must have opposite signs, as can be seen in Fig. 7. When λ L 23 ðM Z Þ is taken in the interval ½−0.2; 0.1, the system remains perturbative up to the grand unification scale, ∼10 15 -10 17 GeV, even when jλ L 33 ðM Z Þj, jλ R 23 ðM Z Þj ≈ 1, as displayed in the first six plots in Fig. 8. However, it deteriorates quickly for larger values of jλ L 33 ðM Z Þj and jλ R 23 ðM Z Þj [ Fig. 8(g)]. When jλ L 23 ðM Z Þj is large, the parameter space is quite limited for jλ L 33 ðM Z Þj, which is in the [0.10, 0.15] band (Fig. 7). In this case as well, the system is well behaved up to jλ L 23 ðM Z Þj, jλ R 23 ðM Z Þj ≈ 1 [e.g., Fig. 8(h)], and the situation declines for the values above in that the perturbativity bound is reached below the unification scale [e.g., Fig. 8(i)].
Note that although we perform the Yukawa RG running based on the SM augmented with a TeV-scale leptoquark S 1 all the way up to the UV, these equations are prone to changes above the intermediate symmetry-breaking scale provided that there is one (as in the examples studied in Sec. III), mainly due to contributions from the running of the scalars whose masses are around the intermediate scale.
However, these effects are expected to be minor due to the corresponding beta-function coefficients not being large enough to significantly change the logarithmic RG running [71]. This is even more likely to be the case especially if this symmetry-breaking scale, for instance the scale in our scenario where the Pati-Salam symmetry is spontaneously broken to the symmetry of the SM, is considerably close to the scale of the SO(10) symmetry-breaking scale (as in the second case in Sec. III, namely model A 2 ), suggesting that the supposed modification in the Yukawa running is indeed not an issue of concern since the slow logarithmic running would most likely not significantly alter the outcome in this small interval. Threshold corrections due to the intermediate scale are also known to be subleading to the one-loop running [75]. unification is realized (and as long as the intermediate scale is high enough to evade the proton-decay constraints for the terms not forbidden by any symmetry in the Lagrangian). The main point of our case that we emphasize is that if a SO(10) theory is indeed the UV completion of the SM, then one may naturally anticipate a TeV-scale leptoquark accompanying the SM Higgs field, and this could define the field content up to very high energies. Beyond that, one has the freedom to choose the high-energy particle content and the corresponding potential terms in the Lagrangian that lead to an appropriate symmetry-breaking sequence in which the Yukawa couplings remain in the perturbative realm above the intermediate symmetry-breaking scale. With this in mind, even the other cases, displayed in Figs. 8(g) and 8(i), that suffer from the perturbativity problem at relatively low energies could arguably get a pass, as long as the high-energy content of the theory is chosen such that the intermediate symmetry breaking occurs before the perturbativity bound is reached and such that the couplings remain in the perturbative realm up to the unification scale.
We comment in passing on the unification-scale implications of our model regarding the fermion mass spectrum. It has been known in the literature that obtaining a realistic Yukawa sector in the SO(10) framework is not trivial. None of the single-and dual-field combinations of the scalar fields 10 H , 120 H , and 126 H yields GUT-scale relations between fermion masses consistent with the SM values [68,75]. On the other hand, it was concluded in Ref. [ The results, obtained by using the equations given in the Appendix and the same values for the input parameters at M Z as in Ref. [75], are given in Table V, some of which are displayed in terms of mass ratios relevant in the SO(10) framework. The unification scale is selected as M U ¼ 2 × 10 16 GeV, which is around the exponential midpoint of the unification scales of our models A 1 and A 2 discussed in Sec. III; small numerical differences depending on the M U value in that range are not relevant to our discussion here. Our values for the fermion masses at M U in the SM case differ from the ones in Ref. [75] by 1-5%, which might be due to a combination of effects coming from the running of the right-handed neutrinos above the intermediate scale and the threshold corrections, both of which are ignored in our estimation. 6 The results for the leptoquark case are given for three benchmark points in the ðλ L 33 ; λ R 23 ; λ L 23 Þ parameter space that are consistent with the perturbativity analysis above. As expected, when all three of the new Yukawa couplings are in the deep perturbative region, the differences from the SM values remain insignificant. When the couplings get larger close to unity, some deviations are observed yet they do not become substantial. Therefore, we conclude that the addition of S 1 to the particle content up to the TeV scale does not lead to significant changes regarding the expected fermion masses at the unification scale and therefore the analysis of Ref. [75], based on using the extrapolated SM values at the unification scale as input data in order to numerically fit them to the parameters of the Yukawa sector of their SO (10) model, can be applied in our case as well.

V. SUMMARY AND DISCUSSIONS
In this paper, we considered the scenario that there is a single scalar leptoquark, S 1 , at the TeV scale, and it is the color-triplet component of a real 10 of SO (10), which also contains a SU(2) doublet which is identified as the SM Higgs. In this scenario, the leptoquark being the only scalar entity other than the SM Higgs is natural; a peculiar mass splitting between the components of 10 does not occur and the leptoquark picks up an electroweak-scale mass together with the SM Higgs, as expected. This is appealing because this leptoquark by itself can potentially explain the B-decay anomalies at the LHC [21,22,25,32].
The SO(10) grand unification framework is an appealing scenario, which has been heavily studied in the literature [62,[64][65][66][67][68][69][70][71][72][73][74][75]. It unifies the three forces in the SM, explains the quantization of electric charge, provides numerous dark matter candidates, accommodates a seesaw mechanism for small neutrino masses, and justifies the remarkable cancellation of anomalies through the anomaly-free nature of the SO(10) gauge group. Moreover, the fermionic content of the SM fits elegantly in 16 F , including a right-handed neutrino for each family. The particles in the SM (except probably the neutrinos) acquire masses up to the electroweak scale through interacting with the electroweak-scale Higgs field, which is generally assigned to a real 10 H . Considering that the leptoquark S 1 is the only other component in 10 H , a TeV-scale S 1 , from this perspective, is consistent with the idea that it is the last piece of the  16 GeV that are selected for the sake of argument as the exponential midpoint of the unification scales of the high-energy scenarios discussed in Sec. III. Clearly, inclusion of S 1 in the low-energy particle spectrum does not lead to significant changes in the fermion mass values at the unification scale compared to the SM predictions, especially in the deep perturbative region in the ðλ L 33 ; λ R 23 ; λ L 23 Þ parameter space.  6 We do not discuss neutrino masses in this paper but have no reason to suspect that the outcome would be different than that in Ref. [75] particularly because their intermediate symmetrybreaking scales at which the generation of neutrino masses through the seesaw mechanism occurs are quite close to the ones in our models A 1 and A 2 . puzzle regarding the particle content up to the electroweak scale (modulo the right-handed neutrinos). Therefore, the possible detection of S 1 , in the absence of any other new particles, at the TeV scale could be interpreted as evidence in favor of SO(10) grand unification.
One obvious issue of concern is proton decay since the leptoquark S 1 possesses the right quantum numbers for it to couple to potentially dangerous diquark operators. On the other hand, the proton stability could possibly be ensured through various mechanisms [22,38,68,[76][77][78]. In this paper, we ensured the proton stability by assuming a discrete symmetry that is imposed in an ad hoc manner below the Pati-Salam breaking scale. It would certainly be more compelling to realize a similar mechanism at the fundamental level although it seems unlikely that it would interfere with the bottom line of this work.
Having a single S 1 leptoquark at the TeV scale as the only new physics remnant from our SO(10) GUT model, we investigated how competent this S 1 is at addressing the R D ðÃÞ anomalies while simultaneously satisfying other relevant constraints from flavor, electroweak and direct LHC searches. We adopted a specific Yukawa coupling texture with only three free (real) nonzero couplings viz. λ L 23 , λ L 33 and λ R 23 . We have found that this minimal consideration can alleviate the potential tension between the R D ðÃÞ -favored region and the R νν K ðÃÞ measurements. The Z → ττ decay constrains λ L 33 whereas the ττ resonance search at the LHC puts complimentary bounds on λ L 23 and λ R 23 . By combining these constraints with all the relevant flavor constraints coming from the latest data on R D ðÃÞ , F L ðD Ã Þ, P τ ðD Ã Þ, R νν K ðÃÞ we have found that a substantial region of the R D ðÃÞ -favored parameter space is allowed. Our multiparameter analysis clearly shows that contrary to the common perception, a single leptoquark solution to the observed R D ðÃÞ anomalies with S 1 is still a viable solution.
Evidently, by introducing new d.o.f. into our framework, one can, a priori, enlarge the allowed parameter region. For example, by considering some of the couplings as complex or by choosing a complex (instead of a real) 10 representation of SO(10), which will introduce an additional S 1 and another complex scalar doublet at the TeV scale, one can relax our obtained bounds. We also pointed out search strategies for S 1 at the LHC using symmetric and asymmetric pair-and single-production channels. Systematic studies of these channels are discussed elsewhere [104].