Charged Lepton Flavor Violating Processes and Scalar Leptoquark Decay Branching Ratios in the Colored Zee-Babu Model

We consider a neutrino mass generating model which employs a scalar leptoquark, $\Delta$, and a scalar diquark, $S$. The new scalars $\Delta$ and $S$ carry the standard model $SU(3)_c\times SU(2)_L\times U(1)_Y$ quantum numbers $(3,1,-1/3)$ and $(6,1,-2/3)$, respectively. The neutrino masses are generated at the two-loop level, as in the Zee-Babu model\cite{Zee-Babu}, and $\Delta/S$ plays the role of the doubly/singly charged scalar in the Zee-Babu model. With a moderate working assumption that the magnitudes of the six Yukawa couplings between $S$ and the down-type quarks are of the same order, strong connections are found between the neutrino masses and the charged lepton flavor violating processes. In particular, we study $Z\rightarrow \overline{l} l'$, and $l\rightarrow l' \gamma$ and find that some portions of the parameter space of this model are within the reach of the planned charged lepton flavor violating experiments. Interesting lower bounds are predicted that $B(Z\rightarrow \overline{l} l')\gtrsim 10^{-16}-10^{-14}(10^{-14})\times(1\mbox{TeV}\cdot m_S/7 m_\Delta^2)^2$ and $B(l\rightarrow l' \gamma)\gtrsim 10^{-17}-10^{-16}(10^{-18}-10^{-16})\times(1 \mbox{TeV}\cdot m_S/7 m_\Delta^2)^2$ for neutrino masses being the normal(inverted) hierarchical pattern. The type of neutrino mass hierarchy could also be determined by measuring the charged lepton flavor violating double ratios. Moreover, definite leptoquark decay branching ratios are predicted when there is no Yukawa interaction between the right-handed fermions and $\Delta$ ( the branching fraction of $\Delta$ to a charged lepton and a quark is 50\%), which could help refine the collider search limit on the scalar leptoquark mass.


Introduction
It is now well established that at least two of the active neutrinos are massive. New physics beyond the Standard Model(SM) is required to give small but nonzero neutrino masses. A straightforward remedy is adding the right-handed neutrino(s) to the SM so that the active neutrinos can acquire Dirac masses after the SM electroweak symmetry breaking as what other charged fermions do. However, additional mass suppression mechanisms or very tiny Yukawa couplings, 10 −12 × Y top , are required to bring down the resulting Dirac neutrino masses to the sub-eV level. Alternatively, Majorana neutrino masses are sought to alleviate the problem of huge hierarchy among the Yukawa couplings in the Dirac neutrino cases. Whatever the UV origin of Majorana neutrino mass is, the key is to generate the dimension-5 Weinberg effective operator [2], (LH) 2 , where L and H are the SM lepton doublet and the Higgs doublet, respectively, at the low energy. The Weinberg operator conserves the baryon number but violates lepton number by two units. Since all the SM interactions at the low energies conserve both baryon and lepton numbers, the new interactions responsible for generating the Majorana neutrino masses must break the lepton number and the relevant new degree(s) of freedom must carry lepton number. If the relevant new fields also carry nonzero baryon number, there are no tree-level contributions to the Weinberg operator leading to a nature loop suppression to bring down the resulting Majorana neutrino masses. Therefore, leptoquark, a boson which carries both lepton number and baryon number, is one of the well-motivated candidates to generate small Majorana neutrino masses without excessive fine tuning. Moreover, since leptoquark participates strong interaction, it would be interesting that the new particles relevant to the neutrino mass generation mechanism could be directly probed at the hadron colliders. However, it is impossible to generate the desired Weinberg operator by using only one leptoquark because the new interaction vertices always come in conjugated pairs. Something else in the loop(s) which carries baryon number must also be utilized to have zero net baryon number and non-vanishing lepton number at the end. The di-quark, a boson which carries 2/3 of baryon number, is one of the candidates to work with leptoquark for generating the Weinberg operator a . Neutrino masses aside, the leptoquark and di-quark are common in many new physics models where the lepton number or the baryon number is not conserved [5][6][7][8][9][10], such as the grand-unified theories, technicolor and composite models. Yet without positive results, leptoquark and di-quark had been eagerly searched for since the 1980's. At the colliders, the leptoquark and di-quark could be produced and studied directly. However, the decay rates strongly depend on the unknown couplings between the leptoquark/di-quark and the SM fermions. Thus, the bounds are usually given with specific assumptions on their couplings to the SM fermions, see [11][12][13][14][15][16][17][18][19]. On the other hand, flavor changing processes could be mediated by the leptoquark or di-quark at the tree-level. Strong constraint can be indirectly derived from the low energy flavor changing experiments [20].
Recently, an interesting application of utilizing the scalar leptoquark and scalar diquark to generate the neutrino masses was discussed by [21]. In [21], one scalar leptoquark, ∆ with SM SU (3) c ×SU (2) L ×U (1) Y quantum number (3, 1, −1/3), and one scalar di-quark, S with SM quantum number (6, 1, −2/3), were augmented to the SM particle content and the neutrino masses can be generated through the two-loop radiative corrections b . This twoloop mechanism is very similar to that in the Zee-Babu model [1] except that S/∆ replaces the role of the doubly/singly charged scalar in Zee-Babu model. From now on this model is referred as the colored Zee-Babu Model(cZBM). In the cZBM, the resulting neutrino mass matrix pattern and the mixing angles are determined by Y L and Y S , the Yukawa couplings between leptoquark and di-quark and the SM fermions, see Eq.(2.1). Again, Y S and Y L are arbitrary and a priori unknown. To proceed, we consider the case that the symmetric Y S 's are democratic and the magnitudes of the six (Y S ) ij , where i, j = 1, 2, 3 are the flavor indices, are of the same order. This could be realized in the extra-dimensional models with the right-handed down-type quark bulk wave functions cluster together in the extra spatial dimension(s), for applying the geometric setup to generate a special 4dimensional Yukawa pattern see for example [23,24]. With this working assumption and the fact that m b ≫ m s ≫ m d , the Y L can be determined with some reasonable requirements which will be discussed later. To accommodate all the neutrino data, the tree-level flavor violating processes will be inevitably mediated by ∆ with the realistic Y L Yukawa couplings. Moreover, the rates of these resulting tree-level and also those flavor violating processes induced at the loop level must comply with the current experimental bounds. In addition to Y L , ∆ also admits Yukawa couplings, Y R , which couple ∆ to the right-handed leptons and quarks, see Eq.(2.1). Since both Y L and Y R contribute to the tree-level flavor violating processes incoherently, Y R = 0 is assumed to minimize those rates. A comprehensive numerical study is performed to search for the realistic configurations. We find that sizable portion of the realistic solutions overlap with the designed sensitivities of the forthcoming lepton flavor violation experiments. Moreover, for the neutrino masses in both the normal hierarchy and inverted hierarchy, there are interesting and definite lower bounds on B(Z → ll ′ ) and B(l → l ′ γ) which could be falsified in the future. Also, the type of neutrino mass hierarchy can be determined if the charged lepton flavor violating double ratios are measured to be within some specific ranges. If Y R = 0, the model has concrete predictions for the scalar leptoquark decay branching ratios for both neutrino mass hierarchies. This will help refine the collider search limit on the scalar leptoquark mass for the β = 1/2 case.
The paper is organized as follows. A more detailed discussion on the model is given in section 2. In section 3, we study the connection between the neutrino masses and Y L , and the tree-level flavor violating processes as well. The loop-induced flavor violating processes are discussed in section 5. The numerical study are dealt with and discussed in section 5. Finally, the conclusions are summarized in section 6.

Model and neutrino mass
As mentioned in the previous section, the SM is extended by adding S and ∆. After rotating the lepton fields into their weak basis and the quarks into their mass basis, the most general gauge invariant Yukawa interaction associated with S and ∆ is (2.1) where i, j are the flavor indices and the SU (3) indices are suppressed. Apparently Y S is symmetric in the flavor space while there is no such constraints on Y L , Y R , and y ∆ . Moreover, the lagrangian admits a gauge invariant triple coupling term: (µ∆ * ∆ * S + h.c.). As shown in Fig. 1, the neutrino masses will receive nonzero contributions through 2-loop quantum corrections if both Y L and Y S present. If one writes the effective Lagrangian for neutrino masses as − 1 2 ν C Li (M ν ) ij ν Lj , the neutrino mass matrix can be calculated to be Note that the two-loop integral is similar to the one in the Zee-Babu model [25,26]. When the down-type quark mass is much lighter than colored scalars, the integral is flavor inde-pendent and it can be simplified to For the later use, it is convenient to write the neutrino mass matrix in a compact form Qualitatively speaking, the resulting neutrino mass is about One sees that, due to the 2-loop suppression, with a typical values Y L , Y S ∼ 0.01 and µ, M ∼ 1 TeV, the sub-eV neutrino mass can be easily achieved without excessively fine tuning. However, the simultaneous presence of Y L/R and y ∆ leads to tree-level proton decay as pointed out in [27]. A very small y ∆ 11 is enough to avoid the rapid proton decay problem. Alternatively, the y ∆ term can be eliminated by imposing some ad hoc symmetry. For example, this term can be turned off without upsetting all other interactions if some Z 2 parities {−, −, +, +, +, −, +} are assigned to {L, l R , Q, u R , d R , ∆, S}, respectively. Hence, we leave the proton decay problem aside and simply set y ∆ = 0 in this study.
The most general renormalizable scalar potential including S and ∆ is where the trace is over the color indices. The details of the scalar potential are not relevant for the later discussion. We note by passing that only the SM Higgs doublet can acquire a nonzero vacuum expectation value, H = v/ √ 2, and being solely responsible for the electroweak symmetry breaking(EWSB). The tree-level masses of ∆ and S are shifted after EWSB with m 2 ∆ → m 2 ∆ + λ 1 v 2 /2 and m 2 S → m 2 S + λ 2 v 2 /2. To proceed, we need m ∆ and m S after EWSB as input. Since S and ∆ participate strong interactions, they are best searched for at the hadron colliders but so far none has been found yet. Depending on their couplings to the SM fields, some lower bounds on m S and m ∆ were obtained from the null result of collider searches. The current lower bounds on m ∆ are summarized in Table.1. For an E 6 -type diquark, CMS study gives m S > 6TeV [19]. These limits are very Table 1. Summary of leptoquark mass lower bound (in unit of GeV) from direct search with 95% CL. The values in parentheses are for β = 0.5, and β = 1 otherwise. The leptoquark decays branching ratios into lq and νq are denoted as β and (1 − β), respectively, and λ is the Yukawa coupling for lq∆. The leptoquark is assumed to decay into leptons within only one specific generation.  [18] sensitive to the assumptions of decay branching fraction as well as the flavor dependant coupling strengthes. Hence, in the following numerically analysis, we take m S = 7 TeV and m ∆ = 1 TeV as the benchmark values c .
The triple ∆∆S coupling generates 1-loop correction ∼ µ 2 16π 2 log(µ 2 /m 2 X ) to m 2 X where X = S, ∆. For these quantum corrections to be perturbative, one needs roughly |µ 2 log(µ 2 /m 2 X )| ≤ 16π 2 m 2 X d . On the other hand, the dimensionful parameters in the same scalar potential are expected to be around the same order. These considerations led to similar estimations and µ = (0.1 − 1) TeV is assumed in this study.
At the tree-level, the decay channels for leptoquark are ∆ → ℓ i u j and ∆ → ν i d j . For di-quark, it decays into d i d j , and ∆∆ if kinematically allowed. Given that m S , m ∆ ≫ m t , all the SM final states can be treated massless and the decay widthes can be calculated to be (2.10) c Since S and ∆ are also charged under SM SU (3) and U (1)Y , their 1-loop contributions alter the SM hV V ′ couplings where V V ′ = {γγ, γZ, gg}. Following the analysis in [28][29][30] also the data form [31][32][33][34], we find that the corrections to the signal strengths are not significant, 0.96 < µγγ,γZ < 1.2, for m∆ ∈ [1,3] TeV and mS ∈ [6,8] TeV. d From the Eq.(2.7), µ also has a weak lower bound |µ| 10 −6 TeV×(M/TeV) 2 if YL and YS are required to be less than 1.0.

Neutrino Masses and the Tree-level Flavor Violation
As discussed before, it is assumed that there is no hierarchy among the Y S 's. Since m b ≫ m s ≫ m d , the matrix ω can be broken into the leading and sub-leading parts and ω = ω (0) + ω (1) , where It is easy to check that the leading neutrino mass matrix M Hence, at least one of the active neutrinos is nearly massless, ∼ (m d /m b ) × max(m ν ), and the scenario of quasi-degenerate neutrinos is disfavored in the cZBM. At leading order, (Y L ) 11,21,31 do not enter M (0) ν at all. Therefor, for either normal hierarchy (NH) or the inverted hierarchy(IH) type of the neutrino masses, the eigenmasses are • NH: and • IH: where ∆m 2 ≡ m 2 3 − (m 2 1 + m 2 2 )/2. Moreover, the absolute values of neutrino mass can be obtained by plugging in the well determined neutrino data [35] listed in Table 2. For NH, m 2 ∼ 0.00868 eV and m 3 ∼ 0.0496 eV, and for IH, m 1 ∼ 0.0483 eV and m 2 ∼ 0.0492 eV. For both cases, the total sum of neutrino masses automatically agrees with the limit that m ν < 0.23 eV at 95% C.L. from the cosmological observation [36]. Once m 1,2,3 are fixed, the neutrino mass matrix can be worked out reversely by The standard parametrization is adopted that where c ij and s ij represent cos θ ij and sin θ ij , respectively. In the case of Majorana neutrinos, α 21 and α 31 are the extra CP phases that cannot be determined from the oscillation experiments. For simplicity, all Y S 's are assumed to be real and the 2 Majorana CP phases Table 2. The global-fit neutrino data with 1σ deviation [35].
will not be discussed in this paper. The leading order neutrino mass matrix has 5(=6-1) independent entries e . With the democratic Y S assumption, the effective Majorana mass for (ββ) 0ν -decay m ee ∼ 0.0018 eV for the NH case. For the IH case, m ee ∼ 0.0479 eV which is within the sensitivity of the planned (0νββ) detectors with ∼ 1 ton of isotope [37]. Furthermore, the lightest neutrino mass is ∼ O(10 −5 eV) for both IH and NH cases. For a given set at all; they are arbitrary at this level and will be determined in the next order perturbation. This approximation largely saves the work of numerical study and lays out the base for higher order perturbations beyond ω (0) .
The most important next to leading contribution to M ν comes from (Y S ) 13 . If one also perturbs (Y S ) 23,33 around (Y S ) solution to (Y L ) i1 for a given set {(Y S ) 13 , δ 23 , δ 33 } are: With only a handful of free parameters, all the leptoquark left-handed Yukawa can be reasonably determined solely by the neutrino data. However, further checks are needed to determine whether the above solution is phenomenologically viable. Next, the tree-level flavor violation will be discussed. The leptonic and quark flavor violating processes will be generated by exchanging S and ∆ at the tree-level, see Fig.2. Since S, ∆ are heavy, they can be integrated out below the EWSB scale. After Fierz transformation, we obtain where a, b are the color indices.
There are way too many new free parameters and rich phenomenology in the most general model. To simplify the discussion and to extract the essential physics, we consider the case that the new physics has minimal tree-level flavor violation(TLFV). Note that the TLFV contributions from different chiral structures always add incoherently. To minimize the total TLFV we need to suppress the TLFV from each chiral structure as much as possible. Let's concentrate on the purely left-handed operators first. Observe that (1) A trivial flavor violation free solution is that with It is obvious that these kind of solutions allow only one nonzero entry of Y L , as can be easily seen by looking at Fig.2(a). It always leads to 2 massless neutrinos which has been excluded by the current neutrino oscillation data.
only one row or one column of Y L can be nonzero f and the resulting neutrino masses have two zeros again. However, only Y L 's are relevant to the neutrino masses. One can set Y R = 0 to minimize the TLFV, and use the 9 remaining Y L 's to accommodate the neutrino masses. Then the lower bound on each flavor violation process can be found since any nonzero Y R will add to it.
It is very easy to build a model with Y R ≪ 1 or Y R = 0 and we supplement with two examples. Example one is to introduce an extra U (1) x with two SM like Higgs doublets, will kill Y R (and also y ∆ ) but still allow the charged fermions to acquire the Dirac masses from their Yukawa couplings with H 1 or H 2 . There are other issues needed to be considered in this setup. For example whether the U (1) x is global or local and whether it is free of anomaly. But these issues do not concern us since they are not relevant to this study and there are well-known model-building machineries available to deal with these problems. The second example is promoting the 4-dimensional model into a higher-dimensional version. If the wave functions of l R and u R in the extra spacial dimension(s) are well separated, like in [24], or have very little overlapping, like in [23], the resulting Y R is negligible. Anyway, here Y R = 0 is taken as a phenomenology assumption which minimizes the TLFV. Some remarks on the case of Y R = 0 will be discussed in next section.
A model independent analysis of the effective four-fermion operators was done by [20]. The 90% C.L. upper limits on the normalized Wilson coefficient ǫ ijkl (it is not the totally anti-symmetric tensor), for each 4-fermi operator are extracted and listed in Table 3. We have (3.10) For the TLFV mediated by S, the last term in Eq.(3.8), it is best constrained by the neutral meson mixings. Following the convention in [38], the corresponding Wilson coefficients and Table 3.
Before ending this section, we recap the assumptions and discussion so far: • Y S 's are assumed to be democratic and there is no outstanding hierarchy among these Yukawa couplings. This leads to one nearly massless active neutrino and |Y S | 9 × 10 −3 × (m S /7TeV) from the constrains of neutral meson mixings.
• The Yukawa couplings Y R are turned off to minimize the TLFV.
• For a given set of {µ, m S , m ∆ , (Y S ) 13,23,33 } and any one of the Y L 's, all the remaining 8 Y L can be iteratively determined from the absolute neutrino masses and the U P M N S matrix.

Charged lepton flavor violating process at one-loop
In this section, we shall study the charged lepton flavor violating (cLFV) processes ℓ → ℓ ′ γ, Z → ℓ ′l and the like which are induced at the 1-loop level with the leptoquark running in the loop, see Fig.3.

ℓ → ℓ ′ γ
The effective Lagrangian responsible for the cLFV process ℓ → ℓ ′ γ [39,40] is parameterized as For m ′ ℓ ≪ m ℓ , the partial decay width is given as A straightforward calculation yields 3) where the index q sums over q = u, c, t and r q ≡ m 2 q /m 2 ∆ . d ll ′ L can be obtained by simply switching Y L ↔ Y R in the above expression for d ll ′ R . The loop functions are and they take the limits F 1 → 1/12 and F 2 → 7/6 + (2 ln x)/3 when x → 0. Unlike at the tree-level, the contributions to the cLFV processes from Y L and Y R entangle with each other at the loop-level. Since m q F 2 (r q ) ≫ m l F 1 (r q ) (for q = c, t), generally speaking, the last term in Eq.(4.3) which involves both Y L and Y R gives the most important contribution to d ll ′ R g . Therefore, it is expected that by setting Y R = 0 to minimize the TLFV will also reduce the 1-loop cLFV processes in general. In the Y R = 0 case, d ll ′ R dominates the cLFV processes because m ℓ ≫ m ℓ ′ . With τ −1 µ ≈ Γ(µ → eν e ν µ ) = 1 192π 3 G 2 F m 5 µ and Γ τ = 0.002265GeV, the branching ratios for ℓ → ℓ ′ γ are and where a γ q = 1 + 4r q (ln r q + 1) + O(r 2 q ), a γ u ∼ a γ c ∼ 1.0 and a γ t ∼ 0.82. Numerically, we have (4.7)

Anomalous magnetic dipole moment
Similar calculation with little modification can be carried over for the flavor diagonal cases. For the charged lepton, the anomalous magnetic dipole moment is (4.8) Assuming that |(Y L ) lq | 2 ∼ O(1) and m ∆ = 1 TeV, one has △a e ∼ 8.0 × 10 −16 , △a τ ∼ 1.0 × 10 −8 , and △a µ ∼ 3.0 × 10 −11 . Unless m ∆ ≪ 1TeV and all 3 (Y L ) µq are sizable and in phase, △a µ in this model is too small to accommodate the observed discrepancy a exp µ −a th µ = (2.39 ± 0.79) × 10 −9 at 1σ C.L. [41]. Moreover, the model predicts a tiny positive △a e which goes against the direction of the observed value that a exp e − a th e = −10.6(8.1) × 10 −13 at 1σ C.L. [42]. Of course, a much larger △a µ is possible to explain to observed discrepancy between the experimental measured value and the theoretical prediction if Y R = 0.

Electric dipole moments
If Y R = 0, the 1-loop charged lepton electric dipole moment(EDM), d ℓ ∼ Nc , could be large. For m ∆ = 1TeV, |Y L | ∼ |Y R | ∼ 0.01, and the CP phase is of order one, the g Barring the cases of fine-tuned cancelations and the hierarchical Yukawa couplings typical electron EDM is around 10 −24 e-cm which is already 4 orders of magnitude larger than the current limit |d e | < 8.7× 10 −29 e-cm [43]. Then, how to suppress the EDM's in this model will be a pressing theoretical issue. A plain solution is setting m ∆ 100TeV to avoid the too large EDMs but the phenomenology at the low energies are strongly suppressed as well.
On the other hand, if Y R = 0 there is no EDM at the 1-loop level. In fact, the first non-zero EDM contribution we can construct begins at the 3-loop level and it involves both V CKM and U P M N S . An order of magnitude estimate gives: If Y L takes a typical value ∼ 0.01, m ∆ = 1TeV, and the combined CP phase is ∼ O(1), this 3-loop electron EDM is expected to be |d e | 10 −37 e-cm. This upper bound is slightly larger than the estimated SM upper bound for d e but way below the sensitivity of any EDM measurement in the foreseeable future. Consequently, d e is a useful handle to test the Y R = 0 assumption in the cZBM: once the electron EDM was measured to be greater than 10 −37 e-cm, either the Y R = 0 assumption with m ∆ ∼ O(TeV) must be abandoned or more new physics is needed to go beyond the cZBM.

µ − e conversion
The µ − e conversion(MEC) will be mediated by the leptoquark at the tree-level as shown in Fig.2. For Y R = 0, the relevant 4-fermi operator is The cLFV photon dipole operator discussed in the previous section will also contribute to MEC with an expected relative magnitude ∼ (α/16π 2 ) 2 comparing to the tree-level one.
Following the analysis of [40], a more quantitative estimate for the MEC rate is: where Z is the atomic number and N is the neutron number for a certain nucleus. The overall factor C conv depends on the form factors of the nuclei and the momentum of the muon. For instance, C conv ( 48 22 T i) = 1.2 × 10 −3 [20]. As can be seen, the LFV photon dipole indeed has much smaller contribution to the MEC than the tree-level one.

µ → 3e
In this model, there are no tree-level contributions to the cLFV µ → 3e decay. The µ → 3e process is dominated by the cLFV photon dipole transition and its rate is much smaller than B(µ → eγ). As pointed out in [39,40], the ratio of B(µ → 3e) to B(µ → eγ) is basically model-independent: Similarly, with replacing the charged lepton masses, the ratios in the rare tau decays are The decay branching ratios τ → µ + ee and τ → µµe + are negligible because they are doubly suppressed by two cLFV transition vertices.

Z →lℓ ′
The same Feynman diagrams in Fig.3 with photon replaced by Z boson lead to cLFV Z →ll ′ decays. Since Z is massive, it can also admit the vector or axial-vector couplings other than the dipole transition couplings as in the l → l ′ γ cases. The most general gauge invariant Z →ll ′ amplitude is parameterized as: 16) where the 4-momentums are labeled as in Fig.3. From the above parametrization, the branching ratio can be easily calculated to be 17) and the experimentally measured value Γ Z = 2.4952 ± 0.0023GeV [35] is used in our study. The 4 dimensionless coefficients c Z R,L , d Z R,L can be obtained through a lengthy but straightforward calculation. The physics is rather simple and can be understood qualitatively. However, the full analytic results are not very illustrating and will not be presented here h . Let's focus on the Y R = 0 case to simplify the physics discussion. First of all, the h The details will be given in other place. )c Z R and its contribution is totally negligible in this process. The above qualitative understandings agree very well with our full calculation. Hence, only the leading contribution from c Z R is kept in the study. It is more useful to express the final result in the numerical form: where a Z u = a Z c ≃ −0.125 − 0.077i = −0.1468e i31.63 • and a Z t = 1. The imaginary part of a Z u,c comes from the pole of light-quark propagators in the loop when the light quarks are going on-shell in the Z decay. Also note that this cLFV decay branching ratio is around 10 −7 if the absolute square in Eq.(4.18) is of order unit. The ballpark estimate is below but close to the current experimental limits [35,44].
The interference between the sub-diagrams with u(c) and t running in the loop makes the relative phases between a Z u,c and a Z t observable. This physical phase leads to CP violation and in general B(Z →lℓ ′ ) = B(Z →l ′ ℓ). Following [45,46], the CP asymmetries are quantified as: In this model, we have numerically (4.20) where the shorthand notation Y ℓ ′ ℓ q ≡ (Y L ) * ℓ ′ q (Y L ) ℓq . Interestingly, due to the sizable CP phase, the CP asymmetries and the cLFV decay branching ratios are of the same order. Also, for the later convenience, we define B Z ℓℓ ′ ≡ B(Z →lℓ ′ ) + B(Z → ℓl ′ ). Before closing this section, we should point out a simple but useful scaling relationship between Y S and Y L in this model. Recall that the neutrino mass is proportional to Y S Y 2 L . Therefore, if Y S is re-scaled by Y S → λ −2 Y S , then Y L must goes like Y L → λY L to keep the neutrino mass unchanged. After such rescaling, B(ℓ → ℓ ′ γ), MEC, B Z ℓℓ ′ and η ℓℓ ′ go like λ 4 while ǫ ll ′ qq ′ , ∆a l , and EDM go like λ 2 due to their amplitude nature. This scaling relationship largely helps reduce the computer time in finding the realistic configurations. Now we have everything needed for the numerical and phenomenological study.
i Since mu, mc ≪ mZ, they can also be treated as massless particles in this process.

Scanning strategy
As discussed in Sec.3, once the set {µ, m S , m ∆ , (Y S ) 13,23,33 } plus any one out of the 9 Y L 's are fixed, all the remaining 8 Y L 's can be iteratively determined from the absolute neutrino masses and the U P M N S matrix. In our numerical search, we take m ∆ = 1TeV and m S = 7TeV as the benchmark. Moreover, for each configuration, µ is randomly produced within [0.1, 1]TeV. For each search, the neutrino mixings sin 2 θ 12,13,23 , and the Dirac phase δ cp are randomly generated within the 1 sigma allowed range from the global fit, Tab.2.
For simplicity the two Majorana phases are set to be zero. Then the U P M N S matrix can be determined via Eq.(3.5). For a given U P M N S , we still need to know the absolute neutrino eigen-masses in order to obtain the neutrino mass matrix, see Eq.(3.4). As has been discussed, we assume the lightest neutrino mass is zero. Depending on the neutrino mass hierarchy, the other 2 absolute neutrino masses can be determined from the given ∆m 2 21 and ∆m 2 . These 2 mass squared differences are also randomly generated within the 1 sigma allowed range from the global fit. Then the absolute neutrino mass matrix M IH ν (M N H ν ) for the inverted(normal) hierarchy is ready for use. Next, |(Y L ) 13 | is randomly generated as a real number between 10 −10 and 1.0. Because of the scaling relationship discussed in the previous section, we fix |(Y S ) 33 | = 0.0097 j without losing any generality. Then, |(Y S ) 13,23 | are generated within [0.1, 10] × 0.0097 and they must obey 0.1 < |(Y S ) 13 /(Y S ) 23 | < 10.0 to be consistent with our working assumption. The signs of (Y L ) 13 and (Y S ) 13,23,33 are also randomly assigned with equal probabilities being positive or negative. With the above mentioned values, (Y L ) With all Y L 's ready, the randomly generated configuration is further checked to see whether it is viable. A configuration will be accepted if it pass all the following criteria: • All |Y L |'s are less than one so that the model can be calculated perturbatively.
• All the TLFV satisfy the current experimental limits listed in Tab.3.
• All the loop-level cLFV processes must comply with the latest experimental limits k summarized in Tab.4.
All these plots have |(Y S ) 33 | fixed at its maximally allowed value 0.0097. From here, other configurations can be obtained by simply scaling down Y S by Y S → λ −2 Y S ( with λ > 1 ). In response, all the LFV processes branching ratios move up as λ 4l . In some plots, the dashed arrows are put in to guide the reader's eyes and show the drifting direction of the branching ratios during the re-scaling. As the Y S is dialed down, all the points of B Z ll ′ go up along the indicated direction until the B(µ → eγ) hits the current experimental limit. Interestingly, we can predict the upper limits on B Z ll ′ , ranging from 10 −11 to 10 −9 , for the Y R = 0 case. We stress that these upper limits are tied with the Y R = 0 assumption; they could be much larger if Y R = 0. This part of parameter space of cZBM could be probed at the planned TeraZ collider where about 10 12 Z bosons will be produced per year with a few ab −1 luminosity [50,51]. Moreover, this particular assumption will be ruled out if any excess was measured in the future experiment. On the other hand, the lower limits on these LFV processes are rather robust and insensitive to the Y R = 0 assumption. Similarly, interesting upper and lower bounds on B(ℓ → ℓ ′ γ) and the CP asymmetries η ℓℓ ′ can be predicted in cZBM, see Tab.5. Note that the upper bounds on all three B(ℓ → ℓ ′ γ) are just below the current experimental limits. Any improvement in these measurements will cut across the interesting parameter space of cZBM. On the other hand, the cZBM with democratic Y S and m ∆ ∼ O(TeV) can be falsified if no such cLFV processes had been detected above the predicted lower bounds in the future experiments.
l During the scaling, one needs to recheck the configuration is still phenomenologically viable.     Note that the double ratio of any pair of cLFV process branching ratios is invariant under the Y S re-scaling and independent of m ∆ . Our numerical also has concrete predictions for Y R = 0 and these double ratios depend on the neutrino mass pattern, see Fig.6. Therefore they could provide an intriguing mean to determine the neutrino mass hierarchy. In particular, the neutrino mass hierarchy can be unambiguously determined if the measured values fell into any of the decisive windows listed in Tab.6. Most of these interesting double ratio windows are plagued by either small cLFV branching ratios or very limited parameter space. However, R 5 ≡ B(τ → µγ)/B(τ → eγ) and R 7 ≡ B Z µτ /B Z µe look quite promising. In the cZBM, if R 5 is measured in the future rare tau decay experiment to be within 0.03 and 30, the neutrino masses are of NH. If R 7 < 1.0 is measured in the future Z-factory, the neutrino masses are of IH in the cZBM. Even in the worst scenario that none of the measured double ratios overlap with these stated windows, one could still tell which neutrino mass hierarchy is more probable by simple statistics and probability theory. For example, if both R 4 and R 5 were measured to be ∼ 10 3 , then the IH is roughly 4 times more probable than the NH in the cZBM. The above discussion clearly demonstrates that the neutrino oscillation experiments and the cLFV measurements are complimentary to one another to better understand the origin of the neutrino masses.

Double Ratio
R 9 < 0.01 R 9 > 3 × 10 4 Table 6. The definitions of the cLFV double ratios and the ranges which can be used to determine neutrino mass hierarchy.

Leptoquark decay branching ratios
First, some shorthand notations are introduced: where the quark flavors are summed over. For Y R = 0, the SU (2) L symmetry ensures This corresponds to the β = 1/2 case that 50% of the leptoquark decays into a neutrino and a down-type quark. Since the neutrino is hard to be tracked in the detector, we focus on the decay channels with a high energy charged lepton as the primordial final state and define The above defined quantity is clearly independent of m ∆ and Y S re-scaling. The decay branching ratios for leptoquark from our numerical study are displayed in Fig.7. It can be clearly seen in Fig.7 that, depending on the neutrino mass hierarchy, there are special patterns in the leptoquark decay branching ratios. Roughly speaking, for the IH case, the leptoquark decays are either (1) B ∆ e ∼ 1.0 or (2) B ∆ µ ∼ 55% and B ∆ τ ∼ 45%. On the other hand, for the NH case, the B ∆ µ and B ∆ τ are concentrated in the region roughly enclosed by 0.7 B ∆ µ + B ∆ Surprisingly, our numerical study has not found any configuration which has either B ∆ µ ∼ 100% or B ∆ τ ∼ 100%. Dictated by the neutrino oscillation data, the model predicts that the leptoquark can NOT decays purely into the 2nd or the 3rd generation charged leptons. These concrete branching ratios could be used to provide the new benchmark leptoquark mass limits with a better motivation.

Conclusion
We have studied the cZBM which exploits a scalar leptoquark ∆(3, 1, −1/3) and a scalar diquark S(6, 1, −2/3) to generate neutrino masses at the 2-loop level. The neutrino mass matrix element M ν ij ( i,j=1,2,3) is proportional to the product of (Y L ) ik m d k (Y † S ) kk ′ m d k ′ (Y T L ) k ′ j , see Eq.(2.2). The Yukawa couplings Y L and Y S are a priori unknown and arbitrary. To proceed, we have adopted a modest working assumption that all six |Y S | are of the same order. Then the Y L can be iteratively determined owing to the fact that m ∆ , m S ≫ m b ≫ m s ≫ m d . Moreover, the mass of the lightest neutrino is of order 10 −5 eV and the model disfavors the case of nearly degenerate neutrinos. The tree-level flavor violating processes will be inevitably mediated by ∆ or S with the realistic Y L which accommodates the neutrino data. Due to the different chiral structures, the contributions to the flavor violating processes from Y L and Y R do not interfere with each other at the tree-level. We have considered the case that Y R = 0 to minimized the tree-level flavor violating processes( and expect the same to happen for the loop induced cLFV). We also have argued that Y R = 0 is actually favored by the fact that there is no electron EDM has been observed yet. A comprehensive numerical study has been performed to look for the realistic Y L and Y S configurations which pass all the known experimental constraints on the flavor violating processes. The viable configurations were collected and have been used to calculate the resulting 1-loop charged lepton flavor violating Z → ll ′ and l → l ′ γ. Some of the realistic configurations could be probed in the forthcoming cLFV experiments. Interesting and robust lower bounds have been found for these cLFV, see Tab.5. Moreover, the neutrino mass hierarchy could be determined if the measured cLFV double ratio(s) is/are in some specific range(s), see Tab.6. For Y R = 0, ∆ has 50% of chance decaying into a charged lepton and an up-type quark. Specific ratios j B(∆ → l i u j ) for each generation charged lepton l i have been predicted in this model. Given the potential link between the neutrino masses generation and ∆, it seems well-motivated using the predicted leptoquark branching ratios as a benchmark scenario for the future scalar leptoquark search limits.