One colorful resolution to the neutrino mass generation, three lepton flavor universality anomalies, and the Cabibbo angle anomaly

We propose a simple model to simultaneously explain four observed flavor anomalies while generating the neutrino mass at the one-loop level. Specifically, we address the measured anomalous magnetic dipole moments of the muon, $\Delta a_\mu$ , and electron, $\Delta a_e$, the observed anomaly of $b\rightarrow s l^+ l^-$ in the $B$-meson decays, and the Cabibbo-angle anomaly. The model consists of four colorful new degrees of freedom: three scalar leptoquarks with the Standard Model quantum numbers $(3,3,-1/3),(3,2,1/6)$, and $(3,1,2/3)$, and one pair of down-quark-like vector fermion in $(3,1,-1/3)$. The baryon number is assumed to be conserved for simplicity. Phenomenologically viable solutions with the minimal number of real parameters can be found to accommodate all the above-mentioned anomalies and produce the approximate, close to $1\sigma$, neutrino oscillation pattern at the same time. From general consideration, the model robustly predicts: (1) neutrino mass is of the normal hierarchy type, and (2) ${\cal M}^\nu_{ee}\lesssim 3\times 10^{-4}\mbox{eV}$. The possible UV origin to explain the flavor pattern of the found viable parameter space is briefly discussed. The parameter space can be well reproduced within a simple split fermion toy model.


I. INTRODUCTION
Despite of the amazing success of the Standard Model (SM) of particle physics, new physics (NP) beyond the SM is ineluctable because compelling evidence showing that at least two active neutrinos are massive [1]. In addition to the nonzero neutrino masses, several recent experimental measurements prominently deviated from the SM predictions are also suggestive to NP. In particular, • The measured anomalous magnetic moment of muon (g − 2) µ [1,2] differs from the most recent SM prediction [3,4] by an amount of ∼ 3.7σ: where the uncertainty is the quadratic combination of the experimental and theoretical ones.
Another new measurement at the J-PARC [6] is also expected to improve the experimental uncertainty in the near future.
Note that △a Cs e and △a µ have opposite signs. New models  have been constructed to accommodate both △a Cs e and △a µ . Moreover, [31][32][33] also attempt to incorporate the neutrino mass generation with the observed △a e,µ .
However, the most recent α determination by using the rubidium atom [34] yields a different outcome △a Rb e ≃ (+4.8 ± 3.0) × 10 −13 , whose sign is different from the one of △a Cs e . These two highly precise values of α differ by a tantalizing ∼ 5σ. More independent investigations or measurements are required to resolve this tension. At this moment, we should consider both cases in this work.
Whether these 2 − 5 σ anomalies will persist is not predictable; the future improvement on the theoretical predictions 1 and the experimental measurements will be the ultimate arbiters. However, at this moment, it is interesting to speculate whether all the above-mentioned anomalies and the neutrino mass can be explained simultaneously. In this paper, we point out one of such resolutions.
With the addition of three scalar leptoquarks, T (3, 3, −1/3), D(3, 2, 1/6), and S(3, 1, 2/3), and a pair of vector fermion, b ′ L,R (3, 1, −1/3), the plethoric new parameters ( mostly the Yukawa couplings) allow one to reconcile all data contemporarily. Parts of the particle content of this model had been employed in the past to accommodate some of the anomalies. However, to our best knowledge, this model as whole is new to comprehensively interpret all the observed deviations from the SM.
The paper is laid out as follows: in Sec.II we spell out the model and the relevant ingredients.
Following that in Sec.III we explain how each anomaly and the neutrino mass generation works in our model. Next we discuss the various phenomenological constraint and provide some model parameter samples in Sec.IV . In Sec.V we discuss some phenomenological consequences, and the UV origin of the flavor pattern of the parameter space. Then comes our conclusion in Sec.VI. Some details of our notation and the low energy effective Hamiltonian are collected in the Appendix.
The one closely related to our T (3, 3, −1/3) is S * Like most models beyond the SM, the complete Lagrangian is lengthy, and not illuminating. In this work, we only focus on the new gauge invariant interactions relevant to addressing the flavor anomalies. For simplicity, we also assume the model Lagrangian respects the global baryon number symmetry, and both T and S carry one third of baryon-number to avoid their possible di-quark couplings. Moreover, we do not consider the possible CP violating signals in this model.
where all the generation indices are suppressed to keep the notation simple and it should be understood that all the Yukawa couplings are matrices. Moreover, the model allows two kinds of tree-level Dirac mass term: 3 To simplify the notation, we use "{, }" and " [, ]" to denote the SU (2)L triplet and singlet constructed from two given SU (2)L doublets, respectively. Also, "⊙" means forming an SU (2)L singlet from two given triplets; see Appendix for the details. With the introduction of b ′ , the mass matrix for down-quark-like fermions after the electroweak SSB becomes: where Y d is the SM down-quark three-by-three Yukawa coupling matrix in the interaction basis.
Note that M d is now a four-by-four matrix. This matrix can be diagonalized by the bi-unitary where (d 1 , d 2 , d 3 , d 4 ) and (d, s, b, b ′ ) stand for the interaction and mass eigenstates, respectively.
The new notation, d 4 , is designated for the interaction basis of the singlet b ′ L,R , and b ′ is recycled to represent the heaviest mass eigenstate of down-type quark. One will see that the mass and interaction eigenstates of b ′ are very close to each other from the later phenomenology study.
Similarly, the SM up-type quarks and the charged leptons can be brought to their mass eigenstates by U u L/R and U e L/R , respectively 5 . Since λ's are unknown in the first place, one can focus on the couplings in the charged fermion mass basis, denoted as λ T,D,S , which are more phenomenologically useful. However, note that the mass diagonalization matrices are in general different for the left-handed (LH) up-and LH down-quark sectors. If we pick the flavor indices of λ T to label the charged lepton and down quark mass states, the up-type quark in the triplet leptoquark coupling will receive an extra factor to compensate the difference between U d L and U u L . Explicitly, with the four-by-three matrixÃ † As will be discussed in below,Ã is the extended CKM rotation matrix,Ṽ , andÃ decouples. Now, all λ T , λ D , and λ S are three-by-four matrices. 5 Note that the SM neutrinos are still massless at the tree-level.
In the interaction basis, only the LH doublets participate in the charged-current(CC) interaction. Thus, the SM W ± interaction for the quark sector is However, the singlet b ′ L mixes with other LH down-type-quarks and change the SM CC interaction. In the mass basis, it becomes where Therefore, the SM three-by-three unitary CKM matrix changes into a three-by-four matrix in our model. When the b ′ L decouples, the coupling matrix V reduces to the SM V CKM . Instead of dealing with a three-by-four matrix, it is helpful to consider an auxiliary unitary four-by-four matrix To quantify the NP effect, one can parameterize the four-by-four unitary matrix U d L by a unitary three-by-three sub-matrix, U d L3 , and three rotations as: where s i (c i ) stands for sin θ i (cos θ i ), and θ i is the mixing angle between d iL and b ′ L . In this work, we assume there is no new CP violation phase beyond the SM CKM phase for simplicity. Now, Eq.(21) can be parameterized as Again, we use d, s, b, b ′ to denote the mass eigenstates with m d ≃ 4.7MeV, m s ≃ 96MeV, m b ≃ 4.18GeV, and M b ′ , the mass of b ′ , unknown. The null result of direct searching for the singlet b ′ at ATLAS sets a limit that M b ′ > 1.22T eV [101] (by assuming only three 2-body decays: b ′ → W t, bZ, bH ), and similar limits have obtained by CMS [102,103]. We take M b ′ = 1.5TeV as a reference in this paper. Moreover, all the direct searches for the scalar leptoquarks at the colliders strongly depend on the assumption of their decay modes. Depending on the working assumptions, the exclusion limits range from ∼ 0.5TeV to ∼ 1.6TeV [1]. Instead of making simple assumptions, it will be more motivated to associate the leptoquark branching ratios to neutrino mass generation [104] or the b-anomalies [105,106]. In this paper, we take m T,D,S ∼ M LQ = 1TeV as the reference point. And the constraint we obtained can be easily scaled for a different M LQ or Since all the new color degrees of freedom are heavier than T eV , it is straightforward to integrate them out and perform the Fierz transformation to get the low energy effective Hamiltonian, see Appendix B.

A. Neutrino mass
Instead of using the bi-lepton SU (2) singlet and a charged scalar without lepton number as first proposed in Ref. [107], we employ two leptoquarks which carry different lepton numbers to break (a) the lepton number and generate the neutrino mass radiactively. We start with a general discussion on the 1-loop neutrino mass generation. If there are two scalars φ 1,2 which interact with fermion F k and the neutrino via a general Yukawa coupling parameterized as where the fermion F can carry arbitrary lepton number and baryon number (L F , B F ). If the two scalars do not mix, φ 1 /φ 2 can be assigned with the lepton-number and baryon number (L F − 1, L B )/(L F + 1, L B ) and the Lagrangian enjoys both the global lepton-number U (1) L and the global baryon-number U (1) B symmetries 6 . Without losing the generality, F k is assumed to be in its mass eigenstate with a mass m k . If φ 1,2 can mix with each other, the lepton number is broken by two units, and the Weinberg operator [113] can be generated radiatively. Let's denote φ h(l) as the heavier(lighter) mass state with mass m h (m l ), and parameterize their mixing as φ 1 = c α φ l + s α φ h and φ 2 = −s α φ l + c α φ h , where s α (c α ) is the shorthand notation for sin α(cos α) and α is the mixing angle. The resultant neutrino mass from Fig.1(a) can be calculated as which is exact and free of divergence. Note that for the diagonal element, the combination in the bracket should be replaced by 2Re(κ ik λ ik ). When the mixing is small, this result can also be approximately calculated in the interaction basis of φ 1 and φ 2 .
In our model, the mass eigenstate F can be the SM down-type quark or the exotic b ′ , and 6 For the discussion of the pure leptonic gauge symmetry U (1)L, see for example [108][109][110][111][112]. plays the role of φ 1 /φ 2 , as depicted in Fig.1(b). Assume the D-T mixing is small, then for i = j, and 2Re(λ T ik λ D ik ) should be used in the bracket for the diagonal elements. To have sub-eV neutrino masses, we need roughly if b ( ′ ) -quark contribution dominates. More comprehensive numerical consideration with other phenomenology will be given in section IV.
The Feynman diagrams in the mass basis for the anomalous magnetic dipole moment of charged lepton in general cases.
We also start with a general discussion on the 1-loop contribution to (g − 2) ℓ by adding a fermion, F , and a charged scalar, h. The F -ℓ-h Yukawa interaction can be parameterized as where both F and ℓ are in their mass eigenstates. Here, we have suppressed the flavor indices but it should be understood that both y R and y L are in general flavor dependent. Then, the resulting 1-loop anomalous magnetic moment depicted in Fig.2(a,b) can be calculated as where Q F is the electric charge of F , and △a F l (△a h l ) is the contribution with the external photon attached to the fermion (scalar) inside the loop. We keep △a F l (△a h l ) in the integral form since the analytic expression of resulting integration is not illuminating at all. The physics is also clear from the above expression that one needs m F ≫ m l also both y l R and y l L nonzero to make △a e and △a µ of opposite sign. For m F ≫ m l , we have Namely, where Similar calculation leads to a l → l ′ γ transition amplitude: where β h = (m h /m F ) 2 , ǫ is the polarization of the photon, and For l = µ and l ′ = e, the above transition amplitude results in the µ → eγ branching ratio [114] Br and it must complies with the current experimental limit, Br(µ → eγ) < 4.2 × 10 −13 [115], or thus can be ignored. where When a ≃ b, the function K takes a limit If factoring out the M F = M b ′ , the dipole transition coefficients in Eq.(37) are given by The current upper bound of Br(µ → eγ) amounts to a stringent limit that the relevant |λ S λ D | 10 −5 . Instead of making the product of Yukawa couplings small, the µ → eγ transition from D-S mixing, Fig.3, can be simply arranged to vanish if muon/electron only couples to b ′ /b or the other way around.
Modulating by the leptoquark masses, numerically we have either Sol-1 : or Sol-2 : can accommodate the observed central values of △a Cs e [△a Rb e ] and △a µ simultaneously with vanishing Br(µ → eγ). However, as will be discussed later, only Sol-2 is viable to simultaneously accommodate the neutrino data.
The b → sµμ transition can be generated by tree-level diagrams mediated by T Instead, to bypass the stringent experimental bounds and the fine-tuning of the parameters, we go for the 1-loop box diagram contribution, as shown in Fig.5, which requires only four nonzero In the usual convention, the transition is described by a low energy effective Hamiltonian with Ignoring the tau mass in the loop, Fig.5(a), the effective Hamiltonian generated by the boxdiagram can be easily calculated as where The function has a limit G(x = 1) = −1/2, and G → −1/x when x ≫ 1.
The second contribution from the box diagram with T ± 1 3 and up-type quark running in the loop yields In arriving the above expression, we have made use of the unitarity of V , namely, and It is clear that the contribution from Fig.5(b) is dominated by λ T µb . And the relevant Wilson coefficients are determined to be For a typical value of β ′ T = (1.0TeV/1.5TeV) 2 , β ′ T G(β ′ T ) = −0.3677. We use the following values, for muon, and for the electron counter part, from the global fit to the b → sl + l − data [63] 8 .
If we take V tb V * ts = −0.03975, then it amounts to Since the combination in the squared bracket is positive, the product λ T τ b (λ T τ s ) * has to be negative. The constraints from B s -B s mixing and b → sγ will be carefully discussed in Sec.IV.

D. Cabibbo-angle anomaly
From the unitarity of V 4 , it is clear that and the Cabibbo-angle anomaly(CAA) is naturally embedded in this model. Moreover, the most commonly discussed unitarity triangle becomes Similarly, this model also predicts that and all the other SM CKM unitary triangles are no more closed in general.
The matrix elements are easy to read. For example, we have The mixing θ i is expected to be small, so a smaller universal to leading order is expected as well. By using the Wolfenstein parameterization and the central values from global fit [1], we have Therefore, to accommodate the deficit of 1st row CKM unitarity (Eq.(6) and Eq. (57)) we have Anomaly Requirement Remark Cabibbo angle anomaly |s 1 + 0.233s 2 | ≃ 0.039 (7) Eq. (57) TABLE IV. The requirement for explaining each mechanism/anomaly. For illustration, we take the following values: V tb V * ts = −0.03975, M LQ = 1.0 TeV, and M b ′ = 1.5 TeV.

IV. CONSTRAINTS AND PARAMETER SPACE
As discussed in the previous section, this model is capable to address neutrino mass generation, △a e,µ , b → sµµ, and the CAA. For readers' convenience, all the requirements are collected and displayed in Table IV.
In this section, we should carefully scrutinize all the existing experimental limits and try to identify the viable model parameter at the end.

A. Low energy 2q2l effective Hamiltonian
In this model, the minimal set (MinS T ) of λ T parameters for addressing all the anomalies and neutrino mass consists of five elements: The following are the consequences of adding other Triplet Yukawa couplings outside the MinS T : (1) At tree-level, λ T ed leads to B + → π + eµ, B 0 →ēτ , τ → eK, and µ-e conversion. Then λ T tree-level λ T eb also leads to µ-e conversion, but the constraint is weak due to the | V ub | 2 suppression. On the other hand, at 1-loop level, it generates the unfavored b → see transition. Also, note that λ T eb = 0 is not helpful for generating M ν ee , which is crucial for the neutrinoless double beta decay. (4) Together with λ T τ s , required for b → sµµ, any nonzero λ T ℓ i d (i = e, µ, τ ) gives rise to K + → π + νν at the tree-level, and thus strongly constrained. Moreover, λ T τ s and λ T τ d generate the K-K mixing via the 1-loop box diagram, and thus stringently limited. (5) Together with MinS T , the presence There is no such effective operator at tree-level. b We update this value by using the new data B(B + → K + µ + τ − ) < 4.5 × 10 −5 [1]. c We update this value by using the new data B(D + → π + µ + µ − ) < 7.3 × 10 −8 [1]. d We obtain the limit by using the top quark decay width, Γt = 1.42GeV, and B(t → qll ′ ) < 1.86 × 10 −5 [117].  hand, the introduction of λ T τ d generates s → dµµ transition via the box-diagram which is severely constraint by the K L → µµ data. So it has to be small too. (7) In general, adding λ T µs will cause conflict with the precision Kaon data.
From the above discussion, adding any λ T ∈ MinS T requires fine tuning the parameters. For simplicity, we set any triplet Yukawa couplings outside the MinS T to zero. However, we still need to scrutinize all the phenomenological constraint on the minimal set of parameters. All the potential a This is the effective operator to address the R(D ( * ) ) anomaly.  [116]. By using the parameter set example of Eq.(109), the model predictions, with the signs kept, are displayed in the last column. We take M T = 1TeV for illustration, and the values in last two columns scale as (M T /1TeV) 2 . Note that the coefficients for su(c)ν µ τ and su(c)ν µ µ are zero at tree-level.
In the table, 0 * stems from choosing V * cb λ T µb + V * cb ′ λ T µb ′ = 0 to retain the µ-e universality in b → clν transition as discussed in the text. detectable effective operators from tree-level contribution of MinS T are listed in Table V and Table   VI. And one has to make sure all the constraints have to be met.
In addition to the limits considered in Ref [116], one needs to take into account the constraint from the lepton universality tests in B decays [83]. In particular, the µ-e universality in the b → cl i ν(i = e, µ) transition has been tested to ≃ 1% level [118]. The MinS T of λ T introduces two operators, (bγ αL c)(ν µ γ αL µ) and (bγ αL c)(ν τ γ αL µ), where the first one interferes with the SM CC interaction while the second one does not. On the other hand, there are no electron counter parts. Therefore, it is required that the modification to the b → cµν j transition rate due to the two new operators is less than ∼ 2%. Their Wilson coefficients, the third and the fourth entities from the end in Table VI, are both proportional to V * cb λ T µb + V * cb ′ λ T µb ′ . For simplicity, we artificially set this combination to zero to make sure the perfect µ-e universality in b → clν at tree level, such that the ratio of λ T µb /λ T µb ′ is fixed as well. However, if more parameter space is wanted, this strict relationship can be relaxed as long as the amount of µ-e universality violation is below the experimental precision.
Finally, due to the QCD corrections, the semi-leptonic effective vector operator for addressing the b → sµµ anomaly gets about ∼ +10% enhancement at low energy [119]. However, all the treelevel 2q2l vectors operators listed in Table V and Table VI, as the constraint, also get roughly the same enhancement factor. Therefore, we do not consider this RGE running factor at this moment.
Next, we move on to consider the tree-level effects from the doublet leptoquark. The non-zero λ D τ b and λ D µb , required for addressing △a e,µ and neutrino data, lead to the following relevant low energy effective Hamiltonian, and its neutrino counter part as well, see Eq.(B3). However, the constraint on these operators are rather weak [116] and can be ignored.

B. SM Z 0 couplings
Because b ′ L,R are charged under U (1) Y hypercharge, they interact with the Z 0 boson. In the interaction basis 9 , the SM Z 0 interaction for the down quark sector is where g SM 3 ) ≃ −0.423, s W = sin θ W , and θ W is the Weinberg angle. If we denote b ′ as d 4 , then the above expression can be neatly written as When rotating into the mass basis, due to the unitarity of the four-by-four U d L,R , it becomes 9 Here we temporarily switch back to earlier notation that b ′ L,R represent the interaction basis.
whereL = (1 − γ 5 )/2,R = (1 + γ 5 )/2, and Using the CP-conserving parametrization introduced in Eq. (22), we have It is clear that, with the presence of b ′ L , the tree-level Flavor-Changing-Neutral-Current (FCNC) in the down sector is inevitable unless at most one of θ 1,2,3 being sizable. For simplicity, we assume one nonvanishing (U d L ) 4d F , where F could be one of d, s, b, and all the others are zero. Let's focus on that specific non-zero flavor diagonal Z-d F -d F coupling. The mixing with b ′ leads to The introduction of b ′ L,R leads to a robust prediction that (g for that down-type quark at the tree-level. Namely, in this model, ) are smaller than the SM prediction. This remind us the long standing puzzle of the bottom-quark forward-backward asymmetry, A b F B , which is 2.3σ below the SM value [1]. However, if we pick θ 3 to be nonzero, then the CAA cannot be addressed, see Eq.(68). Moreover, from our numerical study, only θ 1 = 0 is viable to satisfy all experimental limits. Thus we set θ 2 = θ 3 = 0. From Eq.(68 ), we have and if we take θ 1 to be positive. This predicts g dL = g SM dL + s 2 1 /2 at tree-level, but with negligible effect. On the other hand, one may wonder whether the introduction of λ T τ b , λ T µb and λ D τ b can lead to sizable non-oblique radiactive Z-b-b vertex corrections and address both the A b F B anomaly and R b with the later one agrees with the SM prediction. To address the A F B b anomaly and R b simultaneously, one needs to increase g 2 bR and decrease g 2 bL at the same time. We perform the 1-loop calculation in the M S scheme and the on-shell renormalization, and obtain the UV-finite result: where β Z = (M LQ /m Z ) 2 . Note the diagrams with Z attached to the lepton in the loop have imaginary parts, and this is due to that the lepton pair can go on-shell. Unfortunately, these loop corrections are too small, |δg b L,R | ∼ O(10 −5 ) × |λ T,D | 2 , to be detectable. From the above, we conclude that, barring the tree-level FCNC Z coupling, both A b F B and R b receive no significant modification in this model. Of course, future Z-pole electroweak precision measurements [120][121][122] will remain the ultimate judge. If the A b F B deviation endures, one must go beyond this model. We note by passing that more complicated model constructions are possible to address the A b F B anomaly. For example, this anomaly can be addressed by adding an anomaly-free set of chiral exotic quarks and leptons [123,124], or the vector-like quarks [91,92,125] to the SM.
The Wilson coefficient can be easily calculated to be where the one-forth in the parenthesis is the contribution from T ± 1 3 . Note that this C BB and the SM one are of the same sign, and it increases △M s , the mass difference between B s andB s . But, the central value of the precisely measured △M s = 17.757(21)ps −1 [86] is smaller than the SM one. On the other hand, the SM prediction has relatively large, ∼ 10% [126,127], uncertainties arising from the hadronic matrix elements. If putting aside the hadronic uncertainty, this tension could be alleviated in this model by the extended CKM, V * ts V tb ⇒ V * ts V tb = (V * ts c 2 − V * tb s 2 s 3 )V tb c 3 , which reduces the SM prediction. However, it does not work because we set θ 3 = θ 2 = 0 as discussed in Sec.IV B. Instead, we use the 2σ range to constraint the model parameters. Following Refs. [81,128], the NP contribution can be constrained to be where the factor 0.8 is the RGE running effect from µ LQ ≃ 1TeV to µ b , and C SM BB (µ b ) ≃ 7.2 × 10 −11 GeV −2 is the SM value at the scale µ b . From the above, we obtain so that Eq.(82) can be inside the 2σ confidence interval.
Together with Eq.(56) and assuming that V * cb λ T µb + V * cb ′ λ T µb ′ = 0, the requirement of the treelevel µ-e universality in b → clν, one sees that From this inequality, M T must be around or smaller than TeV for this model parameter to stay in the perturbative region. However, this statement strongly relies on the SM prediction of △M s and the values of C µ 9,10 .

D. B → K ( * ) νν
In this model, the B → K ( * ) νν transition can be mediated by T at tree-level and described by an effective Hamiltonian where and all other Wilson coefficients are zero. Following [129], the normalized branching ratio for B → K ( * ) νν is given by where the SM contribution C ν SM ≃ −6.35 and it is dominated by the Z-penguin. Using V tb V * ts = −0.03975 and the current 90%C.L. limits R ν K < 3.9 and R ν K * < 2.7 given by [130], we obtain a constraint This inequality alone implies |λ T µb (λ T τ s ) * | < 0.0955, which is slightly stronger than the constraint derived from B(B + → K + µ + τ − ) < 4.5 × 10 −5 , see Table V.
Since we also set λ T ed i = 0(d i = d, s, b, b ′ ), there is no µ → eγ transition at 1-loop level by default. Therefore, we only focus on the constraint from τ → µγ and τ → eγ. The rare τ → µγ transition can be induced when both λ T τ b ′ and λ T µb ′ are nonzero. The 1-loop diagram are shown in Fig.6(a).
The dipole τ → µγ amplitude can be parameterized as where k is the photon momentum transfer. If ignoring the muon mass, the partial decay width is given as [114] Γ Since the leptoquark T (D) only couples to the LH(RH) charged leptons, it contributes solely to d τ µ R(L) . If ignoring the charged lepton masses in the loop, the dipole transition coefficient can be easily calculated to be In the above, β ′ is the electric charge of the scalar(fermion) in the loop, and the loop functions, Similarly, the contribution from the diagram with leptoquark D − 2 3 and b ( ′ ) running in the loop yields As discussed in Sec.III B, we set λ D µb ′ = 0 in ( Sol-2) to avoid the dangerous µ → eγ transition 10 . Then, we obtain and due to the accidental cancellation in the squared bracket, it vanishes in the limit of m b ≪ m D in this case. Therefore, only d τ µ R needs to be taken into account. From τ τ = (290.3±0.5)×10 −15 s[1], the branching ratio of this rare process is The current experimental bound, B(τ → µγ) < 4.4 × 10 −8 [131], sets an upper bound On the other hand, if both λ D τ b ′ and λ D eb ′ present, we have for τ → eγ transition. From B(τ → eγ) < 3.3 × 10 −8 [131], it gives a weak bound for β ′ D = (1.0/1.5) 2 .

F. b → sγ
Similar to the previous discussion on τ → µγ, the b → sγ(g) transition can be induced when both λ T τ b and λ T τ s are nonzero, see Fig.6(b,c). Moreover, the fermion masses, m τ or m ντ , can be ignored, which corresponds to the β T ≫ 1 limit. From Eq. (92), it is easy to see that R S (x) → 1/(12x) and R F (x) → 1/(6x) when x ≫ 1. Therefore, by plugging in the electric charges in the loop, the b → sγ transition amplitude can readily read as The first one-half factor comes from the T − 1 3 Yukawa couplings which associate with the (1/ √ 2) normalization, see Eq. (11), and the factor 2 in the last term comes from the anti-tau contribution.
We also need to consider b → sg transition because the RGE running will generate the b → sγ operator at the low energy. Because the gluon can only couple to the leptoquark, the amplitude where T (a) is the SU (3) c generator.

G. Neutrino oscillation data
As discussed before, the vanishing λ T ed i are preferred by phenomenological consideration. It is then followed by a robust prediction that M ee = 0, and the neutrinoless double beta decay mediated by M ee vanishes as well. Therefore, the neutrino mass is predicted to be the normal hierarchical(NH). Later we should discuss the consequence if this vanishing-λ T ed i assumption is relaxed.
A comprehensive numerical fit to the neutrino data is unnecessary to understand the physics, and it is beyond the scope of this paper as well. For simplicity, we assume all the Yukawa couplings are real, and thus all the CP phases in the neutrino mixings vanish. However, this model has no problem to fit the CP violation phases of any values once the requirement of all the Yukawa couplings being real is lifted. Moreover, to adhere to the philosophy of using the least number of real parameters, only two more Yukawa couplings, λ D τ b ′ and λ D τ b , are introduced to fit the neutrino data 11 . Together with MinS T , we make use of eleven Yukawa couplings, and the complete minimal set of parameter is Assuming that M D ≃ M T ≃ M LQ , the neutrino mass matrix takes the form where . Note that the leptoquark Yukawa couplings are tightly entangled with the neutrino mass matrix. For instance, λ T τ b ′ /λ T µb ′ = M ν eτ /M ν eµ is required by this minimal assumption.
For illustration, we consider an approximate neutrino mass matrix 12 with all elements being real. It leads to the following mixing angles and mass squared differences 11 Note that we have not employed λ T τ s in the neutrino data fitting yet. 12 It is just a randomly generated example for illustration. There are infinite ones with the similar structure.
Note all of the above values, except δ CP , are inside the 1σ best fit (with SK atmospheric data) range for the normal ordering given by [135], This example neutrino mass matrix captures the essential features of the current neutrino data.
A better fitting to the neutrino oscillation data, including the phase, by using more(complex) model parameters is expected.
In order to reproduce the neutrino mass matrix, the second solution to △a e,µ , Eq.
and pass all experimental limits we have considered, the last column in Tables V and VI.
However, this is because we want to use the minimal number of model parameters. For instance, if 13 The values of △ae correspond to the 1 σ boundaries of △a Cs e [△a Rb e ], and the C 9(10) is close to the fitted value by only using the theoretically clean modes in [63]. All others are taken to be their central values. the complex Yukawa is allowed, the degrees of freedom are doubled. Lowering M b ′ and increasing µ 1 can both make λ D µb and λ S eb ′ smaller as well. We have no doubt that a better fitting can be achieved in this model by using more (complex) free parameters. But we are content with the demonstration about the ability of this model to accommodate all the observed anomalies and explain the pattern of the observed neutrino data with minimal number of real parameters.

H. 0νββ decay
The mixing between T and D breaks the global lepton number, see Appendix B. Here, we consider whether the neutrinoless double beta decay can be generated beyond the contribution from the neutrino Majorana mass. From the low energy effective Hamiltonian, together with the SM CC interaction, the lepton number violating charged current operators can induce the 0νββ process via the diagram, see Fig.7(a). By order of magnitude estimation, the absolute value of the amplitude strength relative to the usual one mediated by the neutrino Majorana mass, Fig.7(b), is given by In arriving the above value, we take µ 3 = 0.5keV, p ∼ 1 MeV for the typical average momentum transfer in the 0νββ process and assume M ν ee ∼ 0.01eV for comparison. Thus, this tree-level process mediated by T and D can be ignored.

I. A recap
After taking into account all the phenomenological limits, we found the following simple assignment with minimal number of real parameters, is able to accommodate △a Cs e [△a Rb e ], △a µ , the Cabibbo angle, and b → sµµ anomalies simultaneously. Moreover, the resulting neutrino mass pattern is very close to the observed one.

A. Neutrino mass hierarchy and neutrinoless double beta decay
Because we set λ T ed i = 0(i = d, s, b, b ′ ), the neutrino mass element M ν ee vanishes and the neutrino mass is of the NH type. Since we have used λ S eb ′ to explain the observed △a e , adding λ T eb ′ is the minimal extension to yield a non-zero M ν ee . Together with λ T µb and λ T µb ′ , the augmentation of λ T eb ′ leads to an effective Hamiltonian, where Numerically, if we set s 1 = 0.039 and V ub = 0.00361.
The Wilson coefficient is severely constrained, C µe < 9.61 × 10 −5 , from the experimental limit of µ-e conversion rate [116,136]. Namely, λ T eb ′ (λ T µb ′ ) * 0.07 unless the cancellation is arranged in Eq. (115). Moreover, in this model, the ratio of neutrino mass element M ν ee to M ν eµ , see Eq. (105), if we include a non-zero λ T eb ′ to generate the ee-component of M ν , the neutrino mass is still of the NH type, and roughly |M ν ee | 3 × 10 −4 eV. The precision required is beyond the capabilities of the near future experiments [137].

B. Some phenomenological consequences at the colliders
The smoking gun signature of this model will be the discovery of b ′ and the three scalar leptoquarks. Once their quantum numbers are identified, the gauge invariant allowed Yukawa couplings and the mechanisms to address the anomalies discussed in the paper follow automatically. The collider physics of leptoquarks have been extensively studied before, and thus we do not have much to add. The readers interested in this topic are referred to the comprehensive review [138] and the references therein. In the paper, we should concentrate on the flavor physics at around or below the Z pole.
However, it is worthy to point out the nontrivial decay branching ratios of the exotic color states. If we assume the mixings among the leptoquarks are small, their isospin members should be approximately degenerate in mass. Then, the decays are dominated by 2-body decay with two SM fermions in the final states. From the Yukawa couplings shown in Eq.(109), the decay branching ratio of leptoquarks can be easily read. For T − 1 3 , its decay branching ratios are For T 2 3 and T 4 3 , the corresponding decay branching ratios are Finally, we have for D and S leptoquarks.
The dominate decay modes of b ′ are b ′ → LQ + l, and b ′ → W − u i (u i = u, c, t) through the mixing of V u i b ′ . Comparing to M b ′ , the masses of final state particles can be ignored. The width for b ′ → W u i is simply given by and The Yukawa coupling between b ′ and the SM Higgs is through the b ′ -d mixing. So the resulting Yukawa coupling gets double suppression from the small θ 1 and the ratio of m d /v 0 , and so the b ′ → Hd decay can be ignored. Similarly the decays of b ′ → Z 0 d i can be ignored as well. By plugging in the parameters we found, the total decay width of b ′ is Γ b ′ ≃ 134.0GeV, and the branching ratios are for M LQ = 1TeV and M b ′ = 1.5TeV.
We stress that the above decay branching ratios are the result of using the example parameter set given in Eq. (109). The decay branching ratios depend strongly on the model parameters, and the branching ratio pattern varies dramatically from one neutrino mass matrix to another 14 . However, one can see that the decay pattern of the heavy exotic color states in this example solution is very different from the working assumption of 100% b ′ → tW, Zb, Hb used for singlet b ′ and other assumption used for the leptoquark searches at the colliders.
Before closing this section, we want to point out some potentially interesting FCNC top 3-body decays in this framework. From the example solution, we have With an integrated luminosity of 3ab −1 and CM energy at 13TeV, about 2.5 × 10 9 top quark pair events will be produced at the LHC. Therefore, only ∼ 5 events which include at least one top decaying in these 3-body FCNC are expected. However, the 3-body FCNC t → cl i l j branching ratios change if adopting a different solution. During our numerical study, we observe that in some cases there are one or two of them at the O(10 −7 ) level, which might be detectable at the LHC.
See [139] for the prospect of studying these potentially interesting 3-body top decay modes at the LHC.

C. Flavor violating neutral current processes
A few comments on the data fitting are in order: From Table V, one sees that the fit almost saturates, ∼ 80%, of the current limit on decay branching ratio B + → K + µ + τ − . Moreover, the fit is very close to the 2σ limit from B s −B s mixing. The solution seems to be stretched to the limit, and the discovery of lepton flavor violating signals are around the corners. However, this is the trade-off of using minimal number of parameters to reproduce all the central values. If instead aiming for the 1σ values, both can be reduced by half. In addition, the Yukawa couplings are tightly connected with the neutrino mass matrix. During our numerical study, we observe that D + → π + µ + µ − and B + → K + µ + τ − can be far below their experimental upper bounds while τ → µγ close to the current experimental limit if using some different neutrino mass matrix or relaxing the strict relationship that V * cb λ T µb + V * cb ′ λ T µb ′ = 0. Therefore, the model has vast parameter space to accommodate the anomalies with diversified predictions, and we cannot conclusively predict the pattern of the rare process rates at the moment.
However, because we need λ T τ b λ T τ s = 0 to explain the b → sµµ anomaly, the b → sτ τ transition will be always generated at the tree-level. From the example solution, we have if taking V tb V * ts = −0.03975. The additional 1-loop contribution via the box diagram similar to that of b → sµµ can be ignored. This Wilson coefficient is roughly six times larger than the SM prediction that C bsτ τ SM ∼ −4.3 [140][141][142], and push the decay branching ratio to Br(B s → τ + τ − ) ≃ 1.6 × 10 −5 . Although the above value is still two orders below the relevant experimental upper limit, ∼ O(10 −3 ), for B s → τ + τ − at LHCb [143] and B + → K + τ + τ − at BaBar [144], this interesting b → sτ τ transition could be potentially studied at the LHCb and Belle II [145], or at the Z-pole [146].
This is to compare with the standard effective Hamiltonian where In this model, the only CC operator can be generated at tree-level is O V L , thus C V R = 0, C S R,L = 0, and C T = 0. From the global fit with single operator [147], the R(D ( * ) ) anomaly can be well addressed if C V L ≃ 0.08. However, the requirement to retain the µ-e universality in b → clν strictly limits the parameter space. By using the real MinS of parameters, we found the model can render Not only the subset of parameters are required to explain the anomalies, the nearly vanishing entities in the Yukawa matrices play vital roles to bypass the strong flavor-changing experimental constraints. The next question is how to understand the origin of this staggering flavor pattern.
Usually, the flavor pattern is considered within the framework of flavor symmetries. It is highly nontrivial to embed the flavor pattern we found into a flavor symmetry, and it is beyond the scope of this paper. Alternatively, we discuss the possible geometric origin of the flavor pattern in the extra-dimensional theories [158][159][160][161].
A comprehensive fitting, including all the SM fermion masses and mixings, like [162][163][164][165] is also beyond the scope of this paper. Instead, here we only consider how to generate the required flavor pattern of the leptoquark Yukawa coupling shown in Eq. (112). To illustrate, we consider a simple split fermion toy model [161]. We assume all the chiral fermion wave functions are Gaussian locating in a small region, but at different positions, in the fifth dimension. Moreover, all the 5-dim Gaussian distributions are assumed to share a universal width, σ SF . As for the three leptoquarks, we assume they do not have the zero mode such that their first Kaluza-Klein(KK) mode are naturally heavy. More importantly, this setup forbids the leptoquark to develop VEV and break the SU (3) c symmetry. In addition, the wavefunctions of the first leptoquark KK mode are assumed to be slowly varying in the vicinity of the fermion cluster, and can be approximated as constants. On the other hand, the SM Higgs must acquires the zero mode so that it can develop v 0 and breaks the SM electroweak symmetry. Then, the 4-dim effective theory is obtained after integrating out the fifth dimension. The effective λ s ij (s = T, D, S) Yukawa couplings is determined by the overlapping of two chiral fermions' 5-dim wave functions times the product of the scalarspecific 5-dim Yukawa coupling and the scalar's 5-dim wavefunction, denoted as N s , in the vicinity where the fermions locate. The 4-dim Yukawa coupling is given by where z i is the center location in the fifth dimension of the Gaussian wave function of particle−i.
It is clear that only the relative distances matter, so we arbitrarily set z τ L = 0 for the LH tau. Note that different fermion chiralities are involved for different scalar leptoquark Yukawa couplings. For example, λ T ij is determined by the separation between the corresponding LH down-quark(z d jL ) and the LH lepton(z ℓ iL ), but z l iR and z d jL are involved for λ S ij . For simplicity, we assume the mass and interaction eigenstates coincide for the down-type quarks and the SM charged leptons. We found the flavor structure can be excellently reproduced if the chiral fermion locations in the fifth dimension are in this given split fermion location configuration. It is well-known that the Pauli matrices can serve as the generators for the 2-dimensional SU (2) representations. Namely, t (2) i = σ i 2 , (i = 1, 2, 3), and they satisfy the relationship [t (2) i , t j ] = iǫ ijk t (2) k ( i, j, k = 1, 2, 3 ). Any SU (2) doublet R (2) transforms as R (2) → U (2) ( θ)R (2) , where U (2) ( θ) = exp(i θ · t (2) ) and θ = (θ 1 , θ 2 , θ 3 ) is the transformation angle vector. It is popular to represent the triplet by a bi-doublet, and the 3-dimensional representation is less discussed in the literature. Equivalently, one can also construct the SU (2) invariants intuitively using the 3dimensional representation. This is the notation adopted in this paper, and we think it might be useful to collect some facts here for the reader.
In order to construct an SU (2) singlet from any two SU (2) triplets, R One can prove that is an SU (2) singlet.
So far, we have suppressed the flavor indices to avoid the unnecessary notational burden, but it is easy to trace and put them back in when needed. For example, the effective Hamiltonian generated by the pair of triplet Yukawa couplings, λ T id and λ T jD , is given by Dγ αL d ē j γ αL e i + 1 2ν j γ αL ν i