Unified Framework for $B$-Anomalies, Muon $g-2$, and Neutrino Masses

We present a model of radiative neutrino masses which also resolves anomalies reported in $B$-meson decays, $R_{D^{(\star)}}$ and $R_{K^{(\star)}}$, as well as in muon $g-2$ measurement, $\Delta a_\mu$. Neutrino masses arise in the model through loop diagrams involving TeV-scale leptoquark (LQ) scalars $R_2$ and $S_3$. Fits to neutrino oscillation parameters are obtained satisfying all flavor constraints which also explain the anomalies in $R_{D^{(\star)}}$, $R_{K^{(\star)}}$ and $\Delta a_\mu$ within $1\, \sigma$. An isospin-3/2 Higgs quadruplet plays a crucial role in generating neutrino masses; we point out that the doubly-charged scalar contained therein can be produced in the decays of the $S_3$ LQ, which enhances its reach to 1.1 (6.2) TeV at $\sqrt s=14$ TeV high-luminosity LHC ($\sqrt s=100$ TeV FCC-hh). We also present flavor-dependent upper limits on the Yukawa couplings of the LQs to the first two family fermions, arising from non-resonant dilepton ($pp \rightarrow \ell^+ \ell^-$) processes mediated by $t$-channel LQ exchange, which for 1 TeV LQ mass, are found to be in the range $(0.15 - 0.36)$. These limits preclude any explanation of $R_{D^{(\star)}}$ through LQ-mediated $B$-meson decays involving $\nu_e$ or $\nu_\mu$ in the final state. We also find that the same Yukawa couplings responsible for the chirally-enhanced contribution to $\Delta a_\mu$ give rise to new contributions to the SM Higgs decays to muon and tau pairs, with the modifications to the corresponding branching ratios being at (2-6)% level, which could be tested at future hadron colliders, such as HL-LHC and FCC-hh.

1 Introduction Among the many reasons to consider physics beyond the Standard Model (SM), an understanding of the origin of neutrino masses stands out, as neutrino oscillations have been firmly established [1] which require nonzero neutrino masses, in contradiction with the SM. While neutrino masses may be accommodated at tree-level simply by the addition of three SM-singlet right-handed neutrino fields having large Majorana masses via the type-I seesaw mechanism [2][3][4][5][6][7], or by the addition of an SU (2) L -triplet scalar (or fermion) via the type-II [7][8][9] (or type-III [10]) seesaw, there are other interesting scenarios where small neutrino masses arise naturally as quantum corrections [11][12][13][14]. These models of radiative neutrino masses, which we focus on in this paper, are more likely to be accessible for direct experimental tests at colliders. (For recent reviews on radiative neutrino mass models and constraints, see Refs. [15,16].) Here we show that the new particles that are present in these models to induce neutrino masses can also play an important role in explaining certain persistent experimental anomalies, viz. the anomalous magnetic moment of the muon (∆a µ ), and the lepton-flavor-universality violating decays of the B meson (R D ( ) and R K ( ) ). There has been a long-standing discrepancy in the measured value of the anomalous magnetic moment of the muon by the E821 experiment at Brookhaven National Laboratory [17] and the SM theory prediction [18], resulting in a value for ∆a µ ≡ a exp µ − a SM µ = (27.4 ± 7.3) × 10 −10 , which indicates a 3.7 σ discrepancy. The muon g − 2 experiment at Fermilab [19] which is currently in the data accumulation stage, in conjunction with more precise calculations of the dominant hadronic vacuum polarization contribution [20][21][22][23][24][25][26], is expected to settle in the near future whether this discrepancy is indeed due to new physics [27]. Meanwhile, it appears to be productive to envision TeV-scale new physics that can account for the observed anomaly. We shall pursue this line of thought here in the presence of an R 2 (3, 2, 7/6) leptoquark (LQ) scalar (in the notation of Ref. [28], where the numbers in parenthesis denote SU (3) c × SU (2) L × U (1) Y quantum numbers) that also takes part in radiative neutrino mass generation.
We show by explicit construction that a model with R 2 (3, 2, 7/6) and S 3 (3, 3, 1/3) LQs plus ∆(1, 4, 3/2) Higgs field can simultaneously explain the R D ( ) , R K ( ) and ∆a µ anomalies, while being consistent with all low-energy flavor constraints, as well as with the LHC limits. A minimal Yukawa flavor structure is suggested that achieves these, while also providing excellent fits to neutrino oscillation parameters. We have also evaluated constraints from √ s = 13 TeV LHC data on the LQ Yukawa couplings to fermions of the first two families arising from non-resonant pp → + i − j processes mediated by t-channel exchange of LQs. These limits on the couplings are found to be in the range (0.15 − 0.36) for a 1 TeV LQ, which would preclude any solution of R D ( ) with new LQ-mediated decays of the B meson involving ν e or ν µ , an a priori logical possibility. We also show that the ∆ ++ scalar from the ∆(1, 4, 3/2) multiplet, which decays to same-sign dileptons for much of the parameter space, can be probed to masses as large as 1.1 TeV at the high-luminosity (HL) phase of the √ s = 14 TeV LHC with 3000 fb −1 of data, as it can be produced via strong interactions in the decay of S 4/3 3 → (R 2 ) −2/3 + ∆ ++ . The mass reach in this new mode is somewhat better than in the standard Drell-Yan (DY) channel. We also find that the same Yukawa couplings responsible for the chirally-enhanced contribution to ∆a µ give rise to new contributions to the SM Higgs decays to muon and tau pairs, with the modifications to the corresponding branching ratios being at a few percent level with opposite signs, which could be tested at future hadron colliders, such as HL-LHC and FCC-hh.
There have been various attempts to explain radiative neutrino masses and a subset of the anomalies in R D ( ) , R K ( ) and ∆a µ using scalar LQs. For instance, Refs. [45,[51][52][53][54][55][56] address neutrino masses and R K ( ) . Similarly, Refs. [57,58] explain radiative neutrino masses, R D ( ) and R K ( ) , while Ref. [59] explains neutrino masses and lepton g − 2. In some cases such explanations are disconnected from neutrino mass generation, in the sense that removing certain particle from the model would still result in nonzero neutrino masses [60,61]. Our approach here is similar in spirit to Ref. [62], which address all three anomalies, viz., R D ( ) , R K ( ) and ∆a µ , in the context of radiative neutrino masses; but unlike Ref. [62] we do not introduce new vector-like fermions into the model. In the model proposed here there is a close-knit connection between the R D ( ) and R K ( ) anomalies, ∆a µ and neutrino mass. In particular, neutrino mass generation requires all particles that play a role in explaining these anomalies. Removing any new particle from the model would render the neutrino to be massless. For other models of radiative neutrino mass using LQ scalars, see Refs. [63][64][65][66][67][68].
The rest of the paper is organized as follows. In Section 2 we present the basic features of the model, including the Yukawa Lagrangian (cf. Section 2.1), scalar potential (cf. Section 2.2), radiative neutrino mass generation mechanism (cf. Section 2.3) and a desired texture for the Yukawa coupling matrices (cf. Section 2.4) consistent with flavor constraints that can explain the flavor anomalies. In Section 3 we discuss how the LQ scalars present in the model explain the R D ( ) and R K ( ) flavor anomalies. In Section 4 we show how the R 2 LQ explains the ∆a µ anomaly. In this section, we also point out the difficulty in simultaneously explaining the electron g − 2 (cf. Section 4.1), as well as the model predictions for related processes, namely, Higgs decay to lepton pairs (cf. Section 4.2) and muon electric dipole moment (cf. Section 4.3). Section 5 summarizes the low-energy constraints on the LQ couplings and masses. Section 6 analyzes the LHC constraints on the LQs. In Section 7 we present our numerical results for two benchmark fits to the neutrino oscillation data that simultaneously explain R D ( ) , R K ( ) and (g − 2) µ anomalies, while being consistent with all the low-energy and LHC constraints. Section 8 further analyzes the collider phenomenology of the model relevant for the ∆ ++ scalar, and makes testable predictions for HL-LHC and future hadron colliders. Our conclusions are given in Section 9.

The Model
The model proposed here aims to explain the B-physics anomalies R D ( ) and R K ( ) , as well as the muon (g − 2) anomaly ∆a µ , and at the same time induce small neutrino masses as radiative corrections. To this end, we choose the gauge symmetry and the fermionic content of the model to be identical to the SM, while the scalar sector is extended to include three new states, apart from the SM Higgs doublet H: Here the numbers within brackets represent the transformation properties under the SM gauge group SU (3) c × SU (2) L × U (1) Y . The superscripts on various fields denote their respective electric charge Q defined as Q = I 3 + Y , with I 3 being the third-component of SU (2) L -isospin. The R 2 and S 3 LQs are introduced to explain R D ( ) and R K ( ) anomalies respectively. The R 2 LQ also explains ∆a µ through a chirally-enhanced operator it induces, which is proportional to the top quark mass. The SU (2) L -quadruplet ∆ field mixes ω 2/3 from R 2 withρ 2/3 from S 3 (the complex conjugate of ρ −2/3 ), which is needed to generate Majorana neutrino masses radiatively. This multiplet, with its characteristic triply-charged component, was introduced to generate tree-level neutrino masses from dimension (d)-7 effective operators in Ref. [69]; here we use it for radiative mass generation, also via d = 7 operators.

Yukawa Couplings
In addition to the SM Yukawa couplings of the fermions involving the Higgs-doublet field H, the following Yukawa couplings of the R 2 and S 3 LQs are allowed in the model: 2 2) Here we have adopted a notation where all fermion fields are left-handed. Q = (u d) T and ψ = (ν e) T are the SM quark and lepton doublets respectively, {i, j} are SU (2) indices, {a, b} are flavor indices, C is the charge conjugation matrix, ij is the SU (2) Levi-Civita tensor, R 2 = iτ 2 R 2 , and τ α (with α = 1, 2, 3) are the Pauli matrices in the doublet representation of SU (2). The color contraction is unique in each term, which is not shown. It is to be noted that S 3 possesses both leptoquark and diquark couplings, as shown in Eq. (2.2), which would lead to potentially dangerous proton decay operators. Therefore, we set the diquark couplingŷ ab to zero in Eq. (2.2), so that baryon number remains unbroken. This is achieved by assigning baryon number B = −1/3 to S 3 and R 2 , along with B = 1/3 for quarks and −1/3 for anti-quarks, and 0 for leptons and anti-leptons.
We redefine fields to go from the flavor basis (u, d, e) to the mass eigenstates (u 0 , d 0 , e 0 ) for the charged fermions (and similarly for the (u c , d c , e c ) fields) via the following unitary rotations in family space: (2. 3) The Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix V CKM is generated in the process and is given by where P , Q are diagonal phase matrices which are unphysical in the SM, but become physical in non-SM interactions, such as the ones involving the LQs. Note that the unitary rotation on the neutrino fields in Eq. (2.3) is the same as for left-handed lepton fields e, and therefore no Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing in the charged weakcurrent interactions of leptons is induced at this stage. For explaining the anomalies in B-decays and in muon g − 2, there is no need to go to the mass eigenstates of the neutrinos; the distinction between the mass and flavor eigenstates will only affect neutrino oscillation phenomenology. For convenience, we also redefine the Yukawa couplings as follows: Eq. (2.2) can now be written in terms of mass eigenstate fermions (except for neutrinos which are still flavor eigenstates) and the redefined Yukawa couplings as Here we have dropped the superscript 0 in the labeling of mass eigenstates. In the discussions that follow, the quark and lepton fields are to be identified as mass eigenstates. Note that the Yukawa coupling matrices f and y, which respectively appear in the d − e c and d − e couplings, also appear in the u − e c and u − e couplings, along with the generalized CKM matrix V . Any texture adopted for f and y should therefor be consistent with flavor violation in both down-type and up-type quark sectors. The flavor indices i and j in f ij (and similarly for f and y) refer to the quark flavor and the lepton flavor respectively. We shall make use of these interactions in explaining the B-anomalies, ∆a µ and radiative neutrino masses.

Scalar Sector
The most general renormalizable Higgs potential involving H, R 2 , S 3 and ∆ is given by: Here {i, j} are SU (2) L indices, {α, β} are SU (3) c indices, τ a are the Pauli matrices, and T a , T a (with a = 1, 2, 3) are the normalized generators of SU (2) in the triplet and quadruplet representations, respectively. 3 Color-singlet contractions not shown explicitly are to be assumed among two colored fields within the same bracket. For example, the λ RS term has the color contraction (R †α 2 τ a R 2α )(S †β 3 T a S 3β ). Here S 3 ij and ∆ ijk are the completely symmetric rank-2 and rank-3 tensors of SU (2), with their components related to those given in Eq. (2.1) as: 3 This potential differs considerably from the one given in Ref. [46], which is missing many terms.
The presence of the quartic coupling with coefficient λ ∆H 3 in Eq. (2.7) will induce a vacuum expectation value (VEV) for the neutral component of ∆, even when µ 2 ∆ > 0 is chosen. The cubic coupling with coefficient µ would then lead to mixing of ω 2/3 andρ 2/3 components of R 2 and S 3 LQ fields. Such a mixing is required to realize lepton number violation and to generate neutrino masses. We shall be interested in the choice µ 2 H > 0 (which leads to electroweak symmetry breaking), and µ 2 R , µ 2 S > 0 (so that electric charge and color remain unbroken), and µ 2 ∆ > 0 -so that ∆ 0 acquires only an induced VEV. To ensure that this desired vacuum is indeed a local minimum of the potential, we now proceed to derive the masses of all scalars in the model.
We denote the VEVs of H 0 and ∆ 0 fields as While v can be taken to be real by a gauge rotation, v ∆ is complex in general. However, all the complexvalued couplings of the potential, i.e. terms in the last line of Eq. (2.7), can be made real by field redefinitions, which we adopt, and consequently minimization of the potential would make v ∆ real as well.
We obtain the following conditions for the potential to be an extremum around the VEVs of Eq. (2.9), assuming that v = 0: We eliminate µ 2 H and µ 2 ∆ using these two conditions. To derive the scalar mass spectrum, we construct the mass matrices from the bilinear terms resulting from expanding the potential in Eq. (2.7) around the VEVs v and v ∆ .
The 2×2 mass matrix involving the mixing of the charge-2/3 LQs in the basis (ω 2/3 ,ρ 2/3 ) is found to be: where (2.14) The mass eigenstates denoted as X 1,2 are given by X 1 = cos ϕ ω 2/3 + sin ϕρ 2/3 , (2.15) X 2 = − sin ϕ ω 2/3 + cos ϕρ 2/3 , where the mixing angle ϕ is defined as The mass eigenvalues of the charge-2/3 LQ fields are then given as The masses for the remaining LQ components (ω 5/3 , ρ 1/3 , ρ 4/3 ) are obtained as follows: As for the ∆ fields, the masses of the triply and doubly-charged components are given by The singly-charged components of H and ∆ will mix, with a mass matrix given by: (2.24) One combination of (H ± , ∆ ± ) fields is the Goldstone boson (G ± ) eaten up by the W ± gauge boson, while the other combination (δ ± ) is a physical charged Higgs field. These fields are with the mass of δ + given by 4 The neutral CP -even scalars do not mix with the CP -odd scalars, since all couplings and VEVs are real. The mass matrix for the CP -even states in the basis (Re H 0 , Re ∆ 0 ) reads as: The resulting mass eignvalues are given by The corresponding mass eigenstates are given by . (2.32) The field h is to be identified as the SM-like Higgs boson of mass 125 GeV. Similarly, the CP -odd scalar mass matrix, in the basis (Im H 0 , Im ∆ 0 ) is given by We identify the Goldstone mode G 0 eaten up by the Z 0 gauge boson and the physical pseudoscalar Higgs boson A 0 as with the mass of A 0 given by The VEV v ∆ must obey the condition v ∆ v from electroweak T -parameter constraint. In presence of v ∆ , the electroweak ρ parameter deviates from unity at tree-level, with the deviation given by [69] δρ Although there are also loop-induced contributions to δρ, arising from the mass splittings among components of ∆, R 2 , S 3 fields which typically have the opposite sign compared to Eq. (2.36), we assume that there is no precise cancellation between these two types of contributions. A parameter ρ 0 , defined as in the MS scheme, θ W being the weak mixing angle, andρ includes leading radiative corrections from the SM), has a global average ρ 0 = 1.00038 ± 0.00020 [1]. Eq. (2.36) can be compared to this global value, with ρ 0 = 1 in the SM, which sets a limit of |v ∆ | ≤ 1.49 GeV, allowing for 3 σ variation.
In the approximation |v ∆ | |v|, one can solve for v ∆ from Eq. (2.11), to get Substituting this into the masses of the Higgs quadruplet components, we obtain [69] where q i is the (non-negative) electric charge of the component field ∆ i (with i = 1, 2, 3, 4 denoting the four components of ∆ given in Eq. (2.1)). We note that there are two possibilities for mass ordering among these components, depending on the sign of the quartic coupling λ H∆ , with m ∆ +++ being either the heaviest or the lightest member. Phenomenology of these scenarios has been studied extensively in Refs. [69][70][71][72]. By choosing all the bare mass parameters µ 2 X (for X = H, R 2 , S 3 , ∆) in Eq. (2.7) to be positive, and the quartic coupling λ H to be positive, the desired minimum can be shown to be a local minimum, as long as the masses of ∆, R 2 , S 3 are well above v 246 GeV. To verify that this minimum is also the absolute minimum of the potential for some range of parameters, further work has to be done, which is beyond the scope of this paper. Since none of the quartic couplings, except for λ ∆H 3 , plays any crucial role for our analysis, it appears possible to achieve this condition. Similarly, there is enough freedom to choose the quartic couplings so that the potential remains bounded from below. We shall discuss below a set of necessary conditions for the potential to be bounded, which will find application in Sec. 4.2 in the discussion of modified rates for h → + − in the model.

Necessary Conditions for Boundedness of the Potential
While the full set of necessary and sufficient conditions on the quartic couplings of Eq. (2.7) for the Higgs potential to be bounded from below is not easily tractable, certain necessary conditions of phenomenological interest (cf. Sec. 4.2) can be analyzed analytically. We focus on the quartic couplings involving only the H and R 2 fields. With SU (2) L and SU (3) c rotations, these fields can be brought to the form where in R 2 , the color indices run horizontally. Here v, x, y can be taken to be real. The quartic terms V (4) (H, R 2 ) can be then written as The necessary and sufficient conditions for boundedness of this potential can now be derived from the co-positivity of real symmetric matrices [73][74][75]: These conditions should hold for any value of the angle α.
Note that from Eq. (2.45) it follows that if (λ HR − λ HR ) is negative, its magnitude cannot exceed about 1.25, since λ H 0.25 is fixed from the mass of h, m h = 125 GeV, if we demand (somewhat arbitrarily) that none of the quartic couplings should exceed about 3.0 in magnitude from perturbativity considerations. This result will be used in the calculation of the modified Higgs branching ratio h → + − in Sec. 4.2.

Radiative Neutrino Masses
Neutrino masses are zero at the tree-level in the model. However, since lepton number is not conserved, nonzero M ν will be induced as quantum corrections. The leading diagrams generating M ν are shown in Fig. 1, mediated by the charge-2/3 LQs. The Yukawa couplings in Eq. (2.6), together with the ∆ R 2 S 3 trilinear term and the ∆ HHH quartic term in the scalar potential (2.7), guarantee lepton number violation. These interactions result in an effective d = 9 operator that violates lepton number by two units, given by O 1 = (ψQ)(ψu c )(HH)H [16,[76][77][78]. Smallness of neutrino mass can be loosely understood even when the new particles have TeV scale masses, owing to a loop suppression factor and a chiral suppression affecting M ν .
The induced neutrino mass matrix arising from Fig. 1 can be evaluated to be where M u = diag{m u , m c , m t } is the diagonal up-quark mass matrix, and κ 1 , κ 2 are respectively the one-loop and two-loop factors given by The leading contribution to M ν is the one-loop term proportional to κ 1 . In evaluating this loop integral we have ignored the masses of the up-type quarks in relation to the LQ masses. In Eq. (2.50) the parameter ϕ is the ω 2/3 −ρ 2/3 mixing angle given in Eq.
(2.17). Since the effective operator for M ν arising from the one-loop diagram is of the type O d=7 eff = ψψHHH † H, which is of d = 7, one should also consider the lower dimensional d = 5 operator O d=5 eff = ψψHH that can be induced at the two-loop level as shown in Fig. 1. In the approximate expression for κ 2 given in Eq. (2.51), the relevant mass scale is that of the heaviest particle in the loop, denoted here by M , defined as M = max(m X 1 , m X 2 , m ∆ 0 ), with m X 1 ,X 2 being the physical masses of the charge-2/3 LQs (cf. Eq. (2.18)) and m ∆ 0 being the physical masses of the quadruplet (cf. Eq. (2.39)). When m X 1 ,X 2 m ∆ 0 , the ratio , which becomes of order unity for m ∆ 0 < 3 TeV or so. However, as we will see later in Section 7, the R 2 LQ is required to have a mass not larger than about 1 TeV in order for it to explain the R D ( ) anomaly. In this case the two-loop diagram is negligible, and therefore, we only include the one-loop contribution in the neutrino fit described in Section 7.2, although the κ 2 term can be important in a more general setting. The overall factor κ 1 in Eq. (2.49) is a free parameter which needs to be of O(10 −8 ) to get the correct order of magnitude for the neutrino masses. Note that the Yukawa matrix elements f ij and y ij must have at least some entries that are of order one in order to explain the B-decay anomalies. κ 1 ∼ 10 −8 can be achieved by taking either the cubic coupling µ in Eq. (2.7) or the induced VEV v ∆ to be small. Both these choices are technically natural, since if either of these parameters is set to zero, lepton number becomes a good symmetry.
We note that the same operator that leads to neutrino masses in this model also induces an effective ∆-quadruplet coupling to the SM leptons. (Recall that ∆ cannot couple to fermions at the tree level in the model.) This can be seen from partner diagrams of Fig. 1, where the SU (2) L components are chosen differently. Ignoring small SU (2) Lbreaking effects, these couplings would all arise from the same effective operator (ψψH † ∆). Therefore, one can write these couplings as being proportional to M ν . Explicitly, we find that the ∆ ++ coupling to leptons has the Yukawa coupling matrix given by where the 1/ √ 3 is a Clebsch-Gordan factor for the ∆ ++ component of the quadruplet in the expansion of the (ψψH † ∆) operator. Eq. (2.52) will play a crucial role in the collider phenomenology of the quadruplet, as discussed in Section 8.

Yukawa Textures
In order to minimize the number of parameters in our numerical fit to R D , R D , R K , R K , (g − 2) µ , and the neutrino oscillation observables, while satisfying all flavor and LHC constraints, we choose the following economical textures for the Yukawa matrices f , f and y defined as in Eq. (2.6) with the first (second) index corresponding to quark (lepton) flavors: Our motivation for the above textures is as follows: Nonzero (f 32 , f 32 ) can explain the anomalous magnetic moment of the muon via chirally-enhanced top-quark loops. The couplings (f 33 , f 22 , f 23 ) are responsible for R D ( ) , while (y 22 , y 32 ) can explain R K ( ) . Similarly, the coupling f 33 is required to suppress the lepton-flavor-violating (LFV) constraint from chirally-enhanced τ → µγ, while simultaneously explaining (g − 2) µ . The remaining parameters (y 23 (33) , y 31 ) in Eq. (2.54) are needed to satisfy the six neutrino oscillation observables (∆m 2 21 , ∆m 2 31 , sin 2 θ 13 , sin 2 θ 23 , sin 2 θ 12 , δ CP ). For more details, see Section 7. We also note that the zeros in the coupling matrices of Eqs. (2.53)-(2.54) need not be exactly zero; but they need to be sufficiently small so that the flavor changing processes remain under control (cf. Section 5).

B-physics Anomalies
In this section, we present our strategy to reconcile the observed tension between experiment and theory in the lepton flavor universality violating observables in the charged-current decays B → D ( ) ν (with = e, µ, τ ) and the neutral-current decays B → K ( ) + − (with = e, µ) by making use of the R 2 and S 3 LQs.
However, the experimental uncertainty on this measurement is rather large at the moment, and any new physics scenario that explains the R D ( ) anomaly automatically explains the R J/ψ anomaly. Therefore, we will not explicitly discuss R J/ψ in what follows.
In order to confront our model with the experimental data for the charged-current processes, we shall consider LQ contributions to the flavor specific process b → cτ −ν . Thus, only the numerator of Eq. (3.1) is modified by the new LQ interactions. To this end, we consider the general low-energy effective Hamiltonian induced by SM interactions as well as the R 2 and S 3 LQs, which is given by where the first term is the SM contribution, while the remaining terms correspond to new physics contribution, with g V,S,T being the Wilson coefficients defined at the appropriate renormalization scale µ R . As shown in Fig. 2, left panel, the ω 2/3 component of the R 2 LQ mediates the b → cτ −ν semileptonic decay via a tree-level contribution. After integrating out the R 2 field, we obtain the following Wilson coefficients at the matching scale µ R = m R 2 : where = e, µ, τ correspond to the outgoing neutrino flavors ν e , ν µ , ν τ respectively. These Wilson coefficients are then run down in momentum to the B-meson mass scale in the leading logarithm approximation, yielding [97] where β (n f ) 0 = 11 − 2n f /3 is the running coefficient, with n f being the number of quark flavors effective in the relevant momentum regime [98,99]. γ S and γ T are anomalous dimension coefficients given by γ S = −8 and γ T = 8/3. Thus, using α s (m Z ) = 0.118, which yields (using QCD running at four loops) α s (m b ) = 0.2167, α s (m t ) = 0.1074 and α s (m R 2 ) = 0.0868, we obtain the following renormalization factors: 5 We see that the tensorial coupling g T becomes less important at µ R = m b , with its value given by . We also note that we have ignored here the mixing between between the Wilson coefficients g S and g T , which is an excellent approximation, as these off-diagonal terms are much smaller than the diagonal terms [101].
The ρ 1/3 component of the S 3 LQ can also contribute in principle to b → cτν via the Wilson coefficient g V given by However, this contribution cannot accommodate R D ( ) as the relevant Yukawa couplings are highly constrained from flavor physics. Any nonzero y 2 is subject to D 0 −D 0 mixing and must be small (cf. Section 5.5), while LHC limits constrain both y 31 and y 32 (cf. Section 6). Furthermore, the product of the Yukawa couplings y 2 and y 3 is strongly constrained by processes such as B → Kνν. It is also worth mentioning that one can induce Wilson coefficient g V of Eq. (3.15) proportional to y 3 y 33 , in conjunction with CKM mixing. However, for = 3, this contribution has an opposite sign compared to the SM, and therefore would require the new contribution to be twice as large as the SM one, bringing it to the non-perturbative regime. For = 1 or 2, there is no interference with the SM term, which would again require large non-perturbative values from the S 3 contribution. Thus we shall ignore these S 3 -induced contributions to R D ( ) . In Section 7.1, we have shown two best fit values of the Yukawa coupling matrices. For these choices of Yukawa couplings, shown in Eqs. (7.3) and (7.4), we get negligible contribution to g V = −5 × 10 −5 for Fit I and g V = 6 × 10 −6 for Fit II from the S 3 LQ, whereas the allowed 1 σ range to explain R D ( ) is [0.072, 0.11]. Therefore, we will only focus on the R 2 contribution to R D ( ) induced through the Wilson coefficients g S and g T . R D and R D induced through the Wilson coefficients g s and g T at µ R = m b with ν τ in the final state are approximately given by [102] R D R SM where the numerical coefficients arise from the relevant form factors. These expressions are applicable for ν e,µ final states as well, but by setting the Re[g τ S ] and Re[g τ T ] terms in Eqs. (3.16) and (3.17) to zero. This is because the new physics and the SM contributions interfere only when ν = ν τ .
The required values for the Wilson coefficient to get a simultaneous fit for both R D and R D is depicted in Fig. 3. We make use of Flavio package [103] that has NNLO QCD and NLO electroweak corrections coded in it, in generating The intersection between the two bands, highlighted by the purple shaded regions, represents the allowed region that satisfies both anomalies. From this plot, we find that Im(g τ S ) must be nonzero, as first noted in Ref. [104], while Re(g τ S ) should be nearly zero to fit R D ( ) . From Eqs. , which is what forces Re(g τ S ) 0. In the right panel, we set Re(g τ S ) = 0, i.e., we set g τ S (or, equivalently, the f 23 coupling) to be purely imaginary, and switch on the f 22 coupling as well, as is the case with our texture in Eq. (2.53). Again, the 1 σ allowed ranges for R D and R D are shown by the The 1 σ allowed ranges for R D and R D in the complex plane of g τ S with g e,µ S = 0. The purple shaded regions correspond to the allowed region that explains both R D and R D . Right: The 1 σ allowed ranges for R D and R D in the plane of (g τ S , g µ S ) (with g e S = 0). The same result is obtained by replacing g µ S with g e S .
light blue and light coral bands, respectively. The same result is obtained by replacing f 22 with f 21 , i.e., by using g e S instead of g µ S . In our numerical fit to R D ( ) in Section 7, we fix m R 2 (f 22 ) close to its minimum (maximum) allowed value from LHC constraints (discussed in Section 6), and find a neutrino mass fit for f 23 and f 33 such that the g µ,τ S values are within the allowed region for both R D and R D ( ) shown in Fig. 3.
The same effective Hamiltonian (3.9) relevant for R D ( * ) also gives rise to the exclusive decay B c → τ ν. Within our model, the branching ratio for this decay is given by [94,105]: (3.18) Here we have used τ [B c ] = (0.507±0.009) ps, f Bc = 0.43 GeV, and m Bc = 6.2749 GeV. The branching ratio BR(B c → τ ν) has not been measured experimentally. Thus, B c lifetime needs to be compared with theoretical calculations [106][107][108][109][110]. With the benchmark fits shown in Sec. 7, we obtain branching ratio at the level of 12 %, which is consistent with the limit quoted in Refs. [102,105,111,112].

Neutral-current Observables
The relevant lepton flavor universality violation ratios R K and R K are defined as The latest LHCb measurement of R K in the q 2 ∈ [1.1, 6] GeV 2 region (q 2 is the invariant mass of the lepton pair in the decays) is [36] R LHCb (3.23) In addition to these LHCb results, Belle has recently announced new measurements on both R K [114] and R K [115], but these results have comparatively larger uncertainties than the LHCb measurements on R K . The effective Hamiltonian describing the new physics contribution to the neutralcurrent process b → sµ + µ − , in presence of S 3 LQ, is given by and C µµ 10 being the Wilson coefficients. Here we have assumed that the new physics couplings to electrons are negligible. We focus on new physics contributions in the b → sµ + µ − channel, i.e. modifying only the numerator of Eq. (3.19). This is motivated by the fact that an explanation of R K ( ) by modifying the b → sµ + µ − decay provides a better global fit to other observables, as compared to modifying the b → se + e − decay [42]. It is known that both R K and R K can be explained by either a purely vectorial Wilson coefficient C µµ 9 < 0, or a purely left-handed combination, C µµ 9 = −C µµ 10 < 0 [47], with the latter combination performing better in the global analysis due to a ∼ 2 σ tension in the BR(B s → µµ) decay which remains unresolved in the C µµ 9 scenario [42]. In our model, the dominant contribution to b → sµ + µ − comes at tree level via the mediation of the ρ 4/3 component of the S 3 LQ, as shown in Fig. 2, right panel. After integrating out the S 3 field, one can extract the Wilson coefficient for b → sµ − µ + decay as: (3.25) The required best fit values of the Wilson coefficients are C 9 = −C 10 = −0.53, with the 1 σ range being [−0.61, −0.45] [42]. In our numerical fit, y 22 and y 32 are fixed by the neutrino mass fit (up to an overall factor), which is then used to fix m S 3 such that the best-fit value of C 9 = −C 10 is obtained from Eq. (3.25). Note that the R 2 LQ can also give rise to b → s + − transition at tree-level with the corresponding Wilson coefficient given by: (3.26) -18 -There is no acceptable fit to R K ( * ) with C 9 = C 10 . Thus, taking the product of couplings f 2α and f 3α to be zero (or very small), one can suppress R 2 contribution to R K ( ) . On the other hand, a loop-level contribution to b → s + − transition can in principle accommodate R K ( ) , but not simultaneously with R D ( ) , due to the stringent limits from τ → µγ [116]. In our numerical fit, therefore, the R 2 contribution will not play a role in explaining R K ( ) .

Muon Anomalous Magnetic Moment and Related Processes
Virtual corrections due to the LQ states can modify the electromagnetic interactions of charged leptons. The contribution from scalar LQ to anomalous magnetic moments has been extensively studied [117][118][119]. In particular, the ω 5/3 component of the R 2 LQ can explain the muon (or electron) anomalous magnetic moment, as it couples to both lefthanded and right-handed fermions, see Eq. (2.6). The new contribution to the anomalous magnetic moment arising from ω 5/3 LQ can be written as [117,120]: where Q q = 2/3 and Q S = 5/3 are respectively the electric charges of the up-type quark and the LQ propagating inside the loop, as shown in Fig. 4. 6 Here x q = m 2 q /m 2 R 2 and we have ignored terms proportional to m 2 /m 2 R 2 in the loop integral. The loop functions appearing in Eq. (4.1) are: Note that the first term in Eq. (4.1) is the LQ contribution to the anomalous magnetic moment without chiral enhancement, whereas the second term is the chirally-enhanced one, which in our case will be proportional to the top-quark mass.

Difficulty with Explaining ∆a e
A discrepancy has also been reported in the anomalous magnetic moment of the electron, denoted as ∆a e , with a somewhat lower significance of 2.4 σ [121]. The signs of ∆a e and ∆a µ are opposite. We have investigated whether ∆a e can also also explained in our framework, but found that the model does not admit a simultaneous explanation of both anomalies, as Figure 4: Chirally-enhanced contribution from the R 2 LQ to the muon anomalous magnetic moment. introducing couplings of the type f αe would lead to a chirally-enhanced contribution to the decay µ → eγ, which is highly constrained. One can attempt to explain both anomalies by simply avoiding chirally-enhanced i → j γ decays by adopting a redefinition of V f ≡ f in Eq. (2.6). However, one introduces V CKM in the down sector leading to strong constraints arising from processes such as K L → e ± µ ∓ , K L → + − , and K −K mixing. A logical option to explain ∆a e would be to choose the Yukawa coupling f 21 to be of O(1), and rely on the charm-quark loop (proportional to f 21 f 21 ), while being consistent with all the flavor constraints and R D ( * ) . However, it turns out that the required values of the Yukawa couplings in this case have been excluded by the latest LHC dilepton constraints on LQ Yukawa couplings and masses from the non-resonant t-channel process pp → + − . These constraints are discussed later in Section 6, and are summarized in Fig. 8. Therefore, simultaneous explanation of the electron and muon anomalous magnetic moments, together with R D ( ) , is not possible in our setup. Thus, we focus on the parameter space required to explain ∆a µ , but not ∆a e , as the former is the more persistent and significant discrepancy. In particular, we set f αe = f αe = 0 in Eq. (2.53) to avoid any ∆a e contribution for our numerical fits discussed in Section 7.

Modified Higgs Decays to Lepton Pairs
The same R 2 LQ interactions that lead to the chirally-enhanced m t /m µ contribution to the muon g − 2 in Fig. 4 will also induce a loop-level correction to the decay of the SM Higgs boson h → µ + µ − . The Feynman diagrams are shown in Fig. 5. In addition to these diagrams which modify the Yukawa couplings directly, one should also take into account correction to the muon mass arising from the R 2 interactions. The relevant diagram is obtained from Fig. 5 by removing the Higgs boson line. The significance of the LQ diagrams in modifying h → µ + µ − decay has been noted recently in Ref. [122]. We have carried out this calculation independently, and found full agreement with the results of Ref. [122]. It is sufficient to compute the coefficient of the d = 6 operator (ψ µL µ R )H(H † H) which is finite, as any loop correction to the d = 4 operator (ψ µL µ R )H will only renormalize the SM operator. The modification to the branching ratio BR(h → µ + µ − ) is found to be (4.6) The loop function F(x, y) can be expanded to first order in y = m 2 t /m 2 R 2 (so that the coefficient of the d = 6 operator is picked out), and also to the required order in x = m 2 h /m 2 t . Although m 2 h /m 2 t ∼ 1, the actual expansion parameter is some factor k times this ratio, with k ∼ 1/10, leading to a rapidly converging series. The function F(x, y) to third order in m 2 h /m 2 t is found to be For our benchmark fits (see Eqs. (7.3) and (7.4)) with m R 2 = 0.9 TeV, the model predictions for µ µ + µ − as a function of the quartic coupling combination (λ HR −λ HR ) is shown in Fig. 6. These predictions are essentially the same for the two benchmark points, so we present our results for Fit I (cf. Eq. (7.3)) in Fig. 6. The coupling λ HR is responsible for the mass splitting between the ω 2/3 and ω 5/3 components of the R 2 LQ (cf. Eqs. (2.13) and (2.19))), which yields a positive contribution to the electroweak ρ-parameter: where N c = 3 for color-triplets like R 2 . Using the current global-fit result for ρ 0 = 1.00038± 0.00020 [1] (with ρ 0 = 1 in the SM) and allowing for 3 σ uncertainty, we obtain an upper bound on the mass splitting ∆m ≤ 55.9 GeV, which yields a corresponding bound on |λ HR | ≤ 1.66. As discussed in Sec. 2.2.1, a necessary condition for the Higgs potential to be bounded from below (cf. Eq. (2.45)) is that for negative values of (λ HR − λ HR ), its magnitude should be below about 1.25. This assumes (somewhat conservatively) that the magnitudes of all quartic couplings lie below 3.0 to satisfy perturbativity. Using the same constraint, we would then have −1.25 ≤ (λ HR − λ HR ) ≤ 4.66 as the preferred range, which is what we shall choose for our numerical study. Our model prediction for µ µ + µ − is shown in Fig. 6 by the solid blue line. We see that the deviation from the SM predictions in this branching is typically at the (2-6)% level. This is fully consistent with the current LHC measurements: µ ATLAS µ + µ − = 1.2 ± 0.6 [123] and µ CMS µ + µ − = 1.19 +0.41 −0.39 (stat.) +0.17 −0.16 (syst.) [124]. For comparison, we quote in Table I the future collider sensitivities for µ µ + µ − from Ref. [125], and the relevant ones are also shown Collider µ µ + µ − µ τ + τ − HL-LHC [126] 9.2% 3.8% HE-LHC [126] 3.4% 2.2% ILC (1000) [127] 12.4% 1.1% CLIC (3000) [128] 11.6% 1.8% CEPC [129] 17.8% 2.6% FCC-hh [130] 0.82% 0.88% Table I: Expected relative precision of the Higgs signal strengths for future colliders. The numbers shown here are for the kappa-0 scenario of Ref. [125].
in Fig. 6 by the horizontal dotted lines. Thus, our predictions for the modified h → µ + µ − signal strength can be tested at the HL-LHC, HE-LHC, as well as at the FCC-hh colliders. It is also worth pointing out that the Yukawa textures needed to simultaneously explain B-anomalies, muon g − 2, and neutrino mass require the f 33 entry to be nonzero, leading to a new contribution to h → τ + τ − . This is also shown in Fig. 6 [132]. For comparison, we quote in Table I the future collider sensitivities for µ τ + τ − from Ref. [125]. Some of these are also shown in Fig. 6 by the horizontal dot-dashed lines. Thus, our predictions for the modified h → τ + τ − signal strength are potentially detectable at future colliders.
As can be seen from Fig. 6, a characteristic feature of the model in the allowed parameter space accessible to future colliders is that while the shift in the branching ratio of h → µ + µ − is downward compared to the SM, it is upward for the branching ratio of h → τ + τ − .

Muon Electric Dipole Moment
LQ interactions can also lead to electric dipole moments (EDM) of the charged leptons (as well as quarks). Existing limits from electron and muon EDMs would place strong constraints on the imaginary part of the Yukawa couplings of the R 2 LQ [133,134]. These constraints are significant only when the LQ couples to both left-and right-handed charged leptons, as depicted in Fig. 4. The lepton EDM arising from these diagrams is given by [117] |d (4.9) In particular, the constraint arising from electron couplings is quite stringent due to the ACME limit |d e | ≤ 1.1×10 −29 e.cm [135]. However, since our model does not give additional contribution to (g − 2) e , we can simply avoid the electron EDM limit by setting the relevant couplings f αe = f αe = 0 in Eq. (2.53). Furthermore, the muon EDM arising from the CKM phase, and from the phases in the matrices P and Q of Eq. (2.4) when varied in their full range [0, 2π], is found to be at most 3 × 10 −22 e-cm, which is well below the current experimental limit of |d µ | ≤ 1.9 × 10 −19 e-cm [136], but may be potentially measurable in future experiments [19,137,138] with high-intensity muon sources [139].

Low-energy Constraints
This section summarizes the most stringent low-energy flavor constraints that are relevant for our model.

α → β γ
These LFV radiative decays arising from LQ loops set some of the most stringent constraints on the couplings of the LQs to µ and τ . As can be seen from Eq. (2.6), the R 2 LQ has both left-and right-handed couplings to charged leptons via the f and f couplings; thus, it can lead to lepton decays both with and without chiral enhancement. The S 3 LQ on the other hand, only couples to left-handed charged leptons, so it cannot induce α → β γ processes with chiral enhancement.
The decay width for the α → β γ mediated by LQ loops is given by [118,120,140] The amplitudes σ R,L arising from the exchange of R 2 LQ can be written as with the loop functions F i (x q ) defined in Eqs. (4.2)-(4.5). Here we generically denote the masses of both 2/3 and 5/3 components of R 2 as m R 2 , assuming them to be degenerate. Note that the amplitude σ q L can be obtained from σ q R with the substitution f ↔ V f . The last terms in Eqs. (5.2) and (5.3) which are proportional to m q are the chirally-enhanced contributions. Similarly, one can obtain the S 3 LQ contribution by replacing the f couplings in the first term of Eq. (5.2) by y, assigning proper charges for the quark (Q q ) and scalar LQ (Q S ), and dropping the f terms in Eq. (5.2).
In the limit m β → 0, which is a very good approximation for both µ → eγ and τ → γ (with = e, µ), and taking into account the u cT f eω 5/3 , u T (V f )e c ω −5/3 , and d T yeρ 4/3 terms in Eq. (2.6), the full expression for α → β γ in our model can be written as Here we have not included the S 3 contribution from theū c L e L ρ 1/3 term, because it is suppressed compared to the d T L ye L ρ 4/3 contribution because of smaller electric charge, as well as due to a CKM-suppression factor and by a Clebsch factor of 2, as can be seen from Eq. (2.6). Similarly, the ω 2/3 component of the R 2 LQ gives sub-dominant contribution proportional to m 2 b /m 2 R 2 compared to the ω 5/3 component, owing to a GIM-like cancellation [64]; so we have not included it in Eq. (5.4). We have displayed the constraint on the Yukawa coupling f from this process in Table II. Table II: Constraints on the Yukawa couplings as a function of LQ mass from α → β γ decay. Constraints on f couplings are obtained by replacing f with (V f ) for the ω 5/3 LQ. Constraints on the S 3 Yukawa coupling y (V y) arising fromd c L e L ρ 4/3 (ū c L e L ρ 1/3 ) are weaker by a factor of 3/2 (6) in comparison to those shown here for the f couplings, suppressed by smaller electric charge and Clebsch factor of 2, as can be seen from Eq. (2.6).

µ − e Conversion
µ − e conversion in nuclei provides a stringent constraint on the product of the Yukawa couplings in our model. The couplings of the S 3 LQ, in conjunction with CKM rotation, is subject to the LFV process from coherent µ − e conversion in nuclei. The branching ratio for this conversion, normalized to muon capture rate, is given by. [16,64,143]: where Γ N is the muon capture rate of the nucleus, p e and E e are respectively the momentum and energy of the outgoing electron, A, Z, and Z eff are atomic number, mass number and effective atomic number of the nucleus, whereas F p is the nuclear matrix element. The experimental limit from gold nucleus provides the most stringent bound [144] of BR < 7.0 × 10 −13 leading to a constraint on the Yukawa coupling: (5.6)

Z → τ τ Decay
Modifications of Z−boson decays to fermion pairs through one-loop radiative corrections mediated by LQs provide another important constraint on the Yukawa couplings of the LQ fileds in the model. We focus our study on the leptonic Z boson couplings as they are the most precisely determined by experiments [1,145]. Within our model, we require the f 33 coupling to be of O(1) to explain the R D ( ) anomaly. Thus we focus on the Z → τ τ decay which provides a constraint of f 33 . The shift in the coupling of τ R with the Z boson arising through loop corrections involving the R 2 LQ is given by [146] Re[δg τ τ R ] = 3|f 33 | 2 16π 2 , and kept terms only to linear orders in these parameters. Using the experimental results on the effective coupling obtained by the LEP collaboration [145], Re[δg τ τ R ] ≤ 6.2 × 10 −4 , we obtain the 1 σ (2 σ) limit on the Yukawa coupling as for the LQ mass of 900 GeV. Within the context of our model and to find a good fit to R D ( ) , we allow this coupling to be in the 2 σ range. A similar constraint on f 32 can be derived, |f 23 | ≤ 1.7 from Z → µ + µ − decay, which is however much weaker than the constraint that one would obtain from τ → µγ, which requires |f 23 f 33 | ≤ 0.3.

Rare D-meson Decays
Rare meson decays also put important constraints on the model parameters. The relevant decays are D 0 → µ + µ − and D + → π + µ + µ − . 7 For effective Lagrangian for these decays mediated by the R 2 and S 3 LQs is given by (cf. Eq. (2.6)) There is also a contribution from the f Yukawa, but it does not come with V CKM rotation, so we do not need to consider this contribution for our choice of f 1α = 0, while deriving the partial decay width for the decay D 0 → µµ. The decay width for D 0 → µµ proportional to the Yukawa couplings f and y is given by From Eq. (5.10), one can obtain the constraint on f 22 using the experimental limit BR(D 0 → µ + µ − ) < 6.2 × 10 −9 [1]: The semileptonic decay D + → π + µµ is mediated by the same term as shown in Eq. (5.9) and we implement the calculation of Ref. [16] to obtain the following decay rate: where the function F is defined as (5.13) 7 In general, the decays B → Kνν and K + → π + νν would provide more stringent constraint on the LQ Yukawa couplings. However, these bounds are avoided by suppressing the LQ couplings to the first generation quarks.
The numerical value of the function F

2.98
GeV. Using f D = 212 MeV, f π = 130 MeV, g D Dπ = 0.59 and the experimental upper limits on the corresponding branching ratio BR(D + → π + µµ) < 7.3 × 10 −8 , we obtain bounds on the f coupling as Similarly, one can find the constraints on Yukawa coupling y 22 , which is weaker by a factor of √ 2 in comparison to f 22 shown in Eqs. (5.11) and (5.14), owing to a Clebsch factor.

D 0 −D 0 Mixing
Both R 2 and S 3 LQs can give rise to D 0 −D 0 mixing via box diagrams. Explicit calculation of the box diagram involving R 2 LQ gives [147] where f D 212 MeV is the D meson decay constant, and C 1 is the Wilson coefficeint given by Here α is the lepton flavor that runs in the box diagrams, which is summed. The renormalized Wilson coefficients C 1 [148][149][150] and the bag factor B 1 [151], evaluated at µ R = 3 GeV scale, are given by From the experimental value |∆m D | = 0.95 +0.41 −0.44 × 10 10 s −1 [1,152], we obtain the limit The same constraint applies to the f coupling as well. However, in addition to the limit quoted in Eq. (5.18), the Yukawa f is also supplemented by Cabbibo rotation, as seen from Eq. (2.6). Thus, for any nonzero entry in the up-sector f 1α or charm-sector f 2α , a nonzero D 0 −D 0 mixing will be induced by the (V f ) term in Eq. (2.6). Consequently, we get a bound on the individual couplings: Similarly, one can obtain a limit on the individual Yukawa y as well, since a nonzero y 1α (or y 2α ) would result in a box diagram contribution to D 0 −D 0 mixing, owing to the CKM mixing. This has contributions from u − ν term in addition to the u − e term in Eq. (2.6). Thus for any nonzero entry in the up-sector or charm-sector in the Yukawa matrix y, the bound is slightly stronger than that shown in Eq. (5.19): It is worth mentioning that the Yukawa couplings y 3α and f 3α also contribute to D-meson mixing. However, these contributions can be safely ignored in the context of our model as they are strongly suppressed by CKM mixing angles by V cb and V ub .

LHC Constraints on Leptoquarks
At the LHC, the R 2 and S 3 LQs can be pair-produced through gg and qq fusion processes, or can be singly produced in association with charged leptons via s-and t-channel quarkgluon fusion processes. The pair production of the LQs at the LHC is solely dictated by the LQ mass, irrespective of their Yukawa couplings, whereas the single production rate depends on both mass and the Yukawa coupling of the LQ. Therefore, the single-production limits are relevant only for larger Yukawa couplings ∼ O(1) [16,153] to the first and secondgeneration quarks. For the benchmark points studied in Section 7, the Yukawa couplings to the first and second generation quarks are not too large (< 1), hence the collider bounds from single-production are not so significant compared to the limits from QCD-driven LQ pair-production. However, we will show in Section 6.2 that there are stringent limits on the Yukawa couplings of the LQ from the the dilepton processes pp → + i − j .

Pair-production Bounds
Once pair-produced at the LHC, each LQ will decay into a quark and a lepton, and the collider limits on these LQ masses depend on the branching ratios to different decay modes.
To impose the bound on the LQ masses, we use the upper limits on the cross-sections from dedicated searches for pair production of first [154,155], second [155][156][157] and third generation [157][158][159] LQs at the LHC and recast them in the context of our model, following the analysis in Ref. [16]. For this purpose, we first implement our model file in FeynRules package [160] and then analyze the signal cross sections using MadGraph5aMC@NLO [161]. Our results for the R 2 LQ are shown in Fig. 7, where the black, red, green, blue, cyan, purple, orange, gray, and brown solid colored lines respectively represent the current bounds from the je, jµ, bτ , tτ , tν, jν, ce, cµ, and jτ decay mode of the LQ. Here the branching ratio of each decay mode is varied from 0 to 1 individually without specifying the other decay modes, which compensate for the missing branching ratios to add up to one. As expected, the bounds on the first and second-generation LQs are much more stringent, as compared to the third-generation case. We will use this information to our advantage while choosing our benchmark points in Section 7. In particular, for the Yukawa ansatz of Eqs. (2.53), the dominant decay modes of the R 2 LQ are: The branching ratios for these decay modes corresponding to the fits presented in Eqs. (7.3) and (7.4) are shown in Table. III. As we can see, the ω 2/3 component of the R 2 LQ dominantly decays to jν and bτ final states, whereas the ω 5/3 component mostly decays to tτ , and jτ final states. Note that the mass of the ω 2/3 component cannot be very different from that of the ω 5/3 component due to the electroweak precision constraints, and hence,  Figure 7: Summary of the updated direct limits from LQ pair-production searches at the LHC for different quark-lepton decay channels of the R 2 LQ. The branching ratio for a specific decay channel of the LQ as indicated in the figure is varied from 0 to 1, while the other decay channels not specified compensate for the missing branching ratios to add up to one. These limits are independent of the LQ Yukawa coupling.  we consider them to be almost degenerate in our analysis. Given the branching ratios in Table. III, the bbτ + τ − final state gives the most stringent constraint on the R 2 LQ mass, which is required to be larger than 859 GeV, as can be seen from Fig. 7. As for the S 3 LQ relevant for R K ( ) anomaly, it can in principle decay to all quark and lepton flavors, due to the CKM-rotations involved in Eq. (2.6). However, the dominant decay modes of the S 3 LQ corresponding to the Yukawa ansatz in Eqs. (7.3) and (7.4) are ρ 4/3 →sµ + , ρ 1/3 →cµ + ,sν , ρ −2/3 →cν . In addition, for m R 2 , m ∆ < m S 3 , the S 3 LQ can decay to the R 2 LQ and the quadruplet scalar ∆, mediated by the trilinear coupling µ in Eq. (2.7) that is responsible for neutrino mass in our model. For our numerical analysis, we focus on the scenario with the R 2 (S 3 ) LQ mass around ∼ 1 TeV (2 TeV) and the quadruplet mass also around 1 TeV. In this case, the S 3 → R 2 + ∆ decay is the dominant one with ∼ 100% branching ratio. In this case, the various components of S 3 decay as follows:

Model
As a consequence, limits on the S 3 LQ mass from the standard LHC searches are not applicable to our scenario. See Section 8 for more details on the S 3 decay signatures at the LHC. For this decay to occur, S 3 mass should exceed that of R 2 LQ.

Dilepton Bounds
Apart from the direct LHC limits from LQ pair-production, there also exist indirect limits from the cross section measurements on the dilepton process pp → + i − j , which could get significantly modified due to a t−channel LQ exchange for large Yukawa couplings. Ref. [47] had derived indirect limits on the LQ mass and Yukawa couplings involving the τ lepton using the previous resonant dilepton searches at the LHC. Meanwhile, a dedicated search [162] for the non-resonant signals in dielectron and dimuon final states has been performed at the √ s = 13 TeV LHC with integrated luminosity 139 fb −1 , which is more appropriate for the t-channel LQ search. Therefore, we use this recent non-resonant dilepton study to derive new indirect limits on the LQ mass and Yukawa couplings. For this analysis, we first implement our model file in FeynRules package [160], then analyze the cross section for pp → + i − j signal using MadGraph5aMC@NLO [161] and compare the quoted observed limits [162] on the cross-section to derive the limits on the Yukawa coupling for a given LQ mass. Our results are shown in Fig. 8 for different Yukawa couplings f iα and f jα (with i = 1, 2; j = 1, 2, 3; α = 1, 2) of the R 2 LQ. Similar bounds can also be derived for the S 3 LQ. There are no bounds on the f 31 and f 32 couplings quoted in Fig. 8, because they involve top-quark initial states, whereas the bounds on f 31 and f 32 come from bottomquark-initiated processes (cf. Eq. (2.6)). Similarly, we do not report any bounds on the Yukawa couplings involving τ -flavor, as there is no corresponding non-resonant dilepton analysis involving taus available so far. Based on the previous analysis [47], we anyway expect the tau-flavor limits to be weaker than the ones quoted here. Note that the bounds derived in Fig. 8 are independent of the LQ branch ratios, unlike the direct limits shown in Fig. 7. As can be seen from Fig. 8

Numerical Fit
In this section, we present our numerical results for the model parameter space that explains the anomalies in R D ( ) , R K ( ) , and ∆a µ within their 1 σ measured values, while being consistent with all the low-energy and LHC constraints discussed above. It is beyond the scope of this work to explore the entire parameter space of the theory; instead we implement all the constraints and find a few benchmark points to explain the anomalies. First of all, we fix the R 2 LQ mass at 900 GeV to satisfy the LHC bound obtained from pair-produced ω 2/3 decaying to bbτ + τ − (cf. Fig. 7 and Table III). Note that m R 2 needs to be around 1 TeV to explain R D ( ) ; making it larger would require larger f 33

Fit to R D ( )
In Fig. 9, we show the allowed parameter space to explain R D ( ) at 1 σ (orange shaded) and 2 σ (light blue shaded) CL in the most relevant Yukawa coupling plane Im(f 23 ) − |f 33 | for a fixed R 2 LQ mass at 900 GeV. We have also fixed f 22 = 0.29, which is the maximum allowed value from the dilepton constraint (see Fig. 8). Note that a nonzero f 22 is required by the neutrino oscillation fit for the textures we have (see Section 7.2), and a larger f 22 helps widen the R D ( ) region. In our numerical analysis to generate Fig. 9, we have made use of the Flavio package [103]. As already noted in Section 3.1 (cf. Fig. 3), the f 23 coupling The horizontal purple band is from the Z → τ τ constraint. The curved green band and cyan bands respectively represent exclusion from LQ pair production in pp → bbτ τ and pp → jjνν channels at LHC. The vertical yellow band corresponds to the exclusion from LFV decay τ → µγ. The dark purple shaded box represents the 1 σ allowed region for R D ( ) that is consistent with all the constraints in this model.
needs to be complex to get a good fit to R D ( ) . Thus, while doing the minimization to get neutrino oscillation fit, we choose the f 23 coupling purely imaginary, as shown in Fig. 9. The dark purple shaded area highlighted in Fig. 9 represents the allowed region that is consistent with all the constraints in our model. The rest of the colored regions are excluded by various constraints discussed in the previous sections. The horizontal purple band is from Z → τ τ constraint (cf. Eq. 5.8). The green and cyan shaded regions respectively represent LHC exclusion from LQ pair-production in bτ and jν decay modes (cf. Fig. 7). The vertical yellow shaded region corresponds to the exclusion from LFV decay τ → µγ (cf. Table II). In the next subsection, we will choose both f 33 and f 23 values from within the allowed region shown in Fig. 9.

Neutrino Fit
In this section, we explicitly show that the neutrino oscillation data can be explained in our model, while being consistent with the B-anomalies and (g −2) µ , as well as satisfying all the where M ν is the diagonal mass matrix and U PMNS is the 3 × 3 PMNS lepton mixing matrix. We numerically diagonalize Eq. (7.2) by scanning over the input parameters with two different textures as shown in Eqs. (2.53) and (2.54). For ease of finding the fits to oscillation data, we factor out m t into the overall factor and define m 0 = m t κ 1 , where κ 1 is given in Eq. (2.50). Furthermore, we perform constrained minimization in which the neutrino observables are restricted to lie within 3 σ of their experimental measured values, for which we use the recent NuFit5.0 values (with SK atmospheric data included) [163]. Our fit results for the two textures given in Eqs. (2.53) and ( Table IV. It is clear that both fits are in excellent agreement with the observed experimental values. The f 33 entry in the benchmark texture shown above is required for fine-tuning at the level of 7% the τ → µγ amplitude arising from top quark loop with a chiral enhancement (cf. Section 5.1). Note that the input parameter f 23 in both Fit I and Fit IIa is purely complex, which is required to get R D ( ) correct (cf. Fig. 9). Furthermore, the same coupling leads to a significant Dirac CP phase, as can be seen from Table IV, consistent with the recent T2K result [166].
Also shown in Table IV are the fit results for R D , R D ( ) , R K ( ) and (g − 2) µ , all of which are within 1 σ of the experimentally allowed range.

Non-standard Neutrino Interactions
The LQs ω 2/3 from R 2 and ρ −2/3 , ρ 1/3 from S 3 have couplings with neutrinos and quarks (cf. Eq. (2.6)). These couplings can induce charged-current NSI at tree-level [16]. Using the effective dimension-6 operators for NSI introduced in Ref. [167], the effective NSI parameters in our model are given by Any non zero entry in the up-sector f 1α and y 1α , relevant for generating tree-level NSI, does not affect the neutrino oscillation fit, as it is suppressed by the up-quark mass. However, Yukawa couplings to the electron and muon sector f 1α and y 1α (α = 1, 2) are highly constrained by the non-resonant dilepton searches at the LHC. The limit on f 11 and f 12 are 0.19 and 0.16, respectively, for 1 TeV LQ mass (cf. Fig. 8). Also, the limit on y 11 and y 12 are 0.16 and 0.15. Thus ε 11 and ε 22 are sub-percent level, and far beyond the reach of forthcoming neutrino experiments. Furthermore, any nonzero y 1α is in conjunction to g g ρ 4/3 Cabibbo rotation and induces (V y) 2α leading to D 0 −D 0 mixing with a constraint given in Eq. (5.20). As noted in Section 6.2, the LHC limits on the LQ Yukawa couplings in the tau sector are weaker, and in principle, one can allow O(1) Yukawa coupling for f 13 and generate a ε 33 which can be as large as 5.6%. However, we require f 23 to be nonzero and O(1) to explain R D ( ) , and the constraint on the product of Yukawa couplings f 13 f 23 is severe due to the D 0 −D 0 bound, see Eq. (5.18). Thus the induced NSI will again be at a sub-percent level. For simplicity, we choose f 1α = y 1α = 0 for all α = 1, 2, 3 (cf. Eq. (2.54)) in both the numerical fits discussed in Section 7.2.

Collider Implications
This model provides an avenue to test a unified description of B-anomalies, muon anomalous magnetic moment and neutrino masses at the LHC through a new decay channel of the S 3 LQ. The presence of the two scalar LQs R 2 and S 3 and the isospin-3/2 scalar multiplet ∆ (especially its triply-and doubly-charged components) give rise to a rich phenomenology for the LHC. In this section, we analyze the production and decay of the doubly-charged component of the scalar multiplet at the LHC and prospective smoking gun signals correlated with the B-anomalies.
Here we propose a unique production mechanism for the doubly-charged scalars at the LHC via the gluon fusion process, as shown in Fig. 10. In the gluon-gluon fusion process, the S 3 LQ can be pair-produced copiously. Once produced, the various components of the S 3 LQ would decay dominantly to the components of the R 2 LQ and ∆ quadruplet, if kinematically allowed (cf. Eq. (6.3)). Here we will mainly focus on the ρ ∓4/3 → ω ±2/3 ∆ ∓∓ decay channel, as ρ 4/3 and ω 2/3 are respectively the components responsible for the R K ( ) and R D ( ) anomalies in our model. Therefore, the signal shown in Fig. 10 provides a direct test of the R K ( ) and R D ( ) explanations at the high-energy LHC.
Another reason we consider the ∆ ±± production via S 3 decay is that the LQ-induced charged-scalar pair-production rate is not as highly suppressed as the DY rate for higher masses. In addition, there will be an enhancement factor for gluon luminosity compared to the quark luminosity, which becomes even more pronounced at higher center-of-mass energies. This can be seen from Fig. 11, where we compare the doubly-charged scalar pairproduction cross-sections at NLO in the DY mode pp → ∆ ++ ∆ −− and in the new LQ mode pp → ∆ ++ ∆ −− + ω 2/3 ω −2/3 (in Fig. 11, ω 2/3 ω −2/3 is collectively denoted as X) for centerof-mass energies √ s =14, 27 and 100 TeV. Note that for the LQ mode, the cross section only depends on the ρ 4/3 LQ mass; however, to make a direct comparison with the DY mode, we have fixed the ω 2/3 mass at 900 GeV (the preferred value for R ( ) D explanation), and for a given ∆ ±± mass in Fig. 11, have chosen the ρ 4/3 mass such that the ρ ∓4/3 → ω ±2/3 ∆ ∓∓ decay branching ratio is ∼ 50% (with the other 50% going to ω ±5/3 ∆ ∓∓∓ ). From Fig. 11, we infer that the production cross-sections for the doubly-charged scalar in the LQ mode are sizable up to the multi-TeV mass range, and the collider reach in the inclusive mode pp → ∆ ++ ∆ −− + X can be significantly enhanced, compared to the pure DY mode (see Section 8.4 for more details).

Decay of Doubly-Charged Scalars
Now we turn to the decay modes of the quadruplet scalar ∆. The doubly charged scalar ∆ ±± can decay to ± ± via the leptonic coupling given by Eq. (2.52). In addition, being a part of the SU (2) L -quadruplet, the covariant derivative term leads to bosonic decay modes (W ± W ± ) of ∆ ±± . On the other hand, when the mass-splitting between consecutive members of the quadruplet are nonzero, cascade decays also open up. One should note that depending on the quartic coupling λ H∆ , there could be two different hierarchies: (a) when λ H∆ > 0, we have m ∆ ±±± < m ∆ ±± < m ∆ ± < m ∆ 0 and (b) when λ H∆ < 0, we have m ∆ ±±± > m ∆ ±± > m ∆ ± < m ∆ 0 (cf. Eq. (2.52)). Therefore, due to mass-splitting, it can decay in cascades via ∆ ±±± X ∓ or ∆ ± X ± (where X = π, W ) depending on whether ∆m > 0 or ∆m < 0. For simplicity, we consider ∆ ±±± to be the lightest member of the ∆ multiplet throughout our analysis. The partial decay widths for different decay modes of ∆ ±± can be written as [71,72]: where the kinematic functions are given by [72] λ(x, y) = 1 + x 2 + y 2 − 2xy − 2x − 2z , (8.7)  Figure 12: Generic decay phase diagram for ∆ ±± in our model, with m ∆ ±± = 1TeV. The dotted, dot-dashed, dashed and thick solid contours correspond to 99%, 90%, 50% and 10% branching ratios respectively for the leptonic, bosonic or cascade decays, whereas ∆m is the mass splitting between the ∆ ++ and the next lightest scalar component.
If ∆ ±± decay to ∆ ± X ± is allowed, the corresponding partial widths will be the same as in Eqs. (8.3)-(8.5). The different scaling factor due to the Clebsch-Gordon coefficient for the quadruplet scalar is taken into account properly for the partial decay width formulae of the doubly charged Higgs given above. For example, the leptonic decay width given in Eq. (8.1) is suppressed by a factor of 2/3, compared to the type-II seesaw scenario [173,174]. On the other hand, the bosonic and cascade decay modes are enhanced by a factor 3/2 in the quadruplet case compared to the triplet scenario [173][174][175]. In Fig. 12, we show the generic decay phase diagram for ∆ ±± in our model, with m ∆ ±± = 1 TeV. The dotted, dot-dashed, dashed and thick solid contours correspond to 99%, 90%, 50% and 10% branching ratios into the leptonic, bosonic or cascade decay modes. The decay phase diagram clearly depicts that the branching ratio to leptonic decay modes of ∆ ±± decreases with v ∆ , whereas the branching ratio to gauge boson decay mode increases with v ∆ . The cross-over happens at v ∆ = 10 −4 GeV with ∆m ∼ 0, similar to the type-II seesaw case [173,174]. As soon as the mass splitting is set to ≥ 10 GeV, cascade decays open up and start dominating depending on the exact value of v ∆ . Note that the mass splitting |∆m| between any two components of ∆ cannot be larger than ∼ 50 GeV due to stringent constraints from electroweak precision data [72].

Comment on 4-body Decay of ∆
In addition to the two-body decays given in Eqs. (8.1)-(8.6), there will also be four-body decay modes of the doubly-charged scalar via the virtual exchange of R 2 and S 3 LQs proportional to the µ term in Eq. (2.7): ∆ ±± → (ω ±2/3 ) (ρ ±4/3 ) , with each LQ decaying to two fermions. These decays will depend on the same parameters that lead to ∆ ±± → ± ± decays. The phase space for these decays would appear to be comparable to the two-body decays, since the latter has a suppression of a loop factor, 1/(16π 2 ) 2 . We have evaluated these four-body decays of ∆ ++ semi-analytically following the procedure outline in Ref. [176], as well as numerically. The two methods gave very similar results. As an example, for a benchmark values of m ∆ ++ = 800 GeV, m R 2 = 1 TeV, m S 3 = 2 TeV, µ = 246 GeV, v ∆ = 10 −4 GeV, and the values of the Yukawa couplings given in Fit I (cf. Eq. (7.3)), the four-body decay width is 2.3 × 10 −15 GeV, which turns out to be much smaller than that for the dileptonic decay, which is 2 × 10 −9 GeV. As v ∆ is increased, the four-body decay may compete with the dileptonic decay; however, in this case ∆ ++ → W + W + decay would dominate. Consequently, the four-body decay of ∆ ++ can be safely ignored in our discussions.

Signal Sensitivity
We focus on the small v ∆ region which gives same-sign dilepton final states from the ∆ ±± decay, because charged leptons with large transverse momenta can be cleanly identified with good resolution and the charge of the leptons can be identified with fairly good accuracy at hadron colliders. For the benchmark fits given in Section 7.2 with normal hierarchy, the dilepton branching ratios of the ∆ ±± → i j for different flavors are as follows: For simplicity, we focus on the µµ final states and consider the signal pp → ∆ ++ ∆ −− +X → µ + µ + µ − µ − + X to derive the sensitivity at future hadron colliders. The relevant SM background is mainly from the multi-top and multi-gauge boson production [177,178]. However, there are several discriminating characteristics of our signal: (a) the invariant mass distributions for same-sign lepton pair from the ∆ ±± decay would peak at a mass value much higher than the SM Z boson mass; and (b) the outgoing leptons will be more energetic compared to the ones produced in the decay of SM gauge bosons, since these leptons are produced from heavy particle ∆ ±± decay. To derive the signal sensitivity, we first implement our model file in FeynRules package [160], then analyze the cross section for the signal using MadGraph5aMC@NLO [161], simulating the hadronization effects with Pythia8 [179] and detector effects with the Delphes3 package [180]. In order to optimize the signal efficiency over the SM background, we impose the following basic acceptance criteria: p T > 15 GeV for each lepton, pseudorapidity |η | < 2.5 and a veto on any opposite sign dilepton pair invariant mass being close to the Z boson mass |M ( + − ) − m Z | > 15 GeV. In addition, events are selected such that the invariant mass for same-sign muon pair is higher than 500 GeV. After passing through all these acceptance criteria, we estimate   Table V: Comparison of the doubly-charged scalar mass reach in the LQ and DY modes (with same-sign di-muon pair final states only) for 3 ab −1 integrated luminosity.
the required luminosities to observe at least 25 events at different center-of-mass energies ( √ s=14, 27, 100 TeV). Our results are shown in Fig. 13. It is clear that for a given luminosity and a given √ s, the doubly-charged scalar mass reach in the LQ mode is higher than that in the DY mode. The mass reach for 3 ab −1 integrated luminosity is summarized in Table V for different center-of-mass energies.
Once we identify the doubly-charged scalar from the multi-lepton signal, the next step is to distinguish the underlying model. In order to identify whether the ∆ ±± 's come from the S 3 LQ decay, accompanied by the ω 2/3 LQs, we can consider the decay chain given in Fig. 10, i.e.
pp → ρ 4/3 ρ −4/3 → ω −2/3 ∆ ++ ω 2/3 ∆ −− → + + − − + τ + τ − + bb . (8.11) In this case, the right combination of the bτ invariant mass peaks at the ω 2/3 LQ mass, if it is produced on-shell from the ∆ decay. Considering the fact that the benchmark fits in our model give 54% branching ratio of ω 2/3 to bτ (cf . Table III), and taking into account the b-tagging and τ -identification efficiencies of ∼ 70% each, we find that at least 25 signal events in the channel given by Eq. (8.11) can be obtained with 3 ab −1 luminosity for the S 3 LQ masses up to 1.5, 2.4 and 5.5 TeV respectively at √ s = 14, 27 and 100 TeV. Hence, it is possible to independently test the unified description of B-anomalies, muon g − 2 and neutrino masses in our model at future colliders.

Conclusion
We have presented a radiative neutrino mass model involving TeV-scale scalar leptoquarks R 2 and S 3 , which can simultaneously explain the R D ( ) , R K ( ) , as well as muon g − 2 anomalies, all within 1 σ CL, while being consistent with neutrino oscillation data, as well as all flavor and LHC constraints. The R 2 LQ is responsible for the R D ( ) and (g − 2) µ , while the S 3 LQ explains the R K ( ) anomaly. The model also features a scalar quadruplet ∆, which is required for the radiative neutrino mass generation. The same trilinear ∆ R 2 S 3 coupling that is responsible for neutrino mass also leads to interesting collider signatures in the S 3 and ∆ decays that can be probed in the forthcoming run of the LHC. Similarly, the same Yukawa couplings responsible for the chirally-enhanced contribution to ∆a µ give rise to new contributions to the SM Higgs decays to muon and tau pairs, with the modifications to the corresponding branching ratios being at 2-6% level, which could be tested at future hadron colliders, such as HL-LHC and FCC-hh.