A unified explanation of $b \to s \mu^+ \mu^-$ anomalies, neutrino masses and $B\rightarrow \pi K$ puzzle

Anomalies in semi-leptonic $B$ decays could indicate new physics beyond the standard model(SM). There is an older puzzle in non-leptonic $B \to \pi K$ decays. The new particles, leptoquarks and diquarks, required to solve the semi-leptonic and the non-leptonic puzzles can also generate neutrino masses and mixing at loop level. We show that a consistent framework to explain the $B$ anomalies and the neutrino masses is possible and we make predictions for certain rare nonleptonic $B$ decays.


I. INTRODUCTION
Searching for beyond the SM (BSM) physics has been the primary focus of the high energy community. Rare B decays have been widely studied to look for BSM effects. Because these decays get small SM contributions, new physics (NP) can compete with the SM and produce deviations from SM predictions. Over the last few years measurements in certain B decays have shown deviations from the SM. These deviations are observed in two groups− in charged current (CC) processes mediated by the b → cτ −ν tansitions and in the neutral current (NC) processes mediated by b → sℓ + ℓ − transition with ℓ = µ, e. We will focus here on the NC anomalies although it is possible that the CC and the NC anomalies are related [1] but we will not explore that possibility here.
While the discrepancies in b → sµ + µ − can be understood with lepton universal new physics [15], hints of LUV in R K and R * K require NP that couple differently to the lepton generations. A well-studied scenario is to assume NP coupling dominantly to the muons though NP coupling to electrons is not ruled out [16,17]. The b → sµ + µ − transitions are defined via an effective Hamiltonian with vector and axial vector operators: dominated nonleptonic b decays and have been studied extensively. The decays in the set include B + → π + K 0 (designated as +0 ), B + → π 0 K + (0+), B 0 → π − K + (−+) and B 0 → π 0 K 0 (00). Their amplitudes are not independent, but obey a quadrilateral isospin relation: Using these decays, nine observables have been measured: the four branching ratios, the four direct CP asymmetries A CP , and the mixing-induced indirect CP asymmetry S CP in B 0 → π 0 K 0 . Shortly after these measurements were first made (in the early 2000s), it was noted that there was an inconsistency among them. This was referred to as the "B → πK puzzle" [37][38][39][40].
Recently the fits were updated [41][42][43]. In Ref. [41] it was observed that the key input to understanding the data was the ratio of the color-suppressed tree amplitude (C ′ ) to the color-allowed (T ′ ) amplitude. Theoretically, this ratio is predicted to be 0.15 < ∼ |C ′ /T ′ | < ∼ 0.5 [44] with a default value of around 0.2. It was found that for a large |C ′ /T ′ | = 0.5, the SM can explain the data satisfactorily. However, with a small, |C ′ /T ′ | = 0.2, the fit to the data has a p value of 4%, which is poor.
needed. The precise statement of the situation is then, the measurements of B → πK decays allow for NP and so in this paper we will assume there is NP in these decays. There are two types of NP mediators that one can consider for the B → πK decays. One is a Z ′ boson that has a flavor-changing coupling tosb and also couples toūu and/ordd.
The second option is a diquark that has db and ds couplings or ub and us couplings. We will focus on the diquark explanation as the diquarks can contribute to neutrino masses.
The paper is organized in the following manner. In Sec. II we describe the setup with leptoquarks and diquarks that leads to neutrino masses and mixing at the loop level. In that section we also discuss the low energy constraints for the leptoquark Yukawa couplings including the b → sℓ + ℓ − data. In Sec. III we explore the B → πK decays mediated by the exchange of diquarks and we consider the constraints on the diquark Yukawa couplings from the B → πK decays and meson oscillations. In Sec. IV we consider the collider constraints on the diquark and leptoquarks coupling and masses and we give a scan of all their couplings that satisfy all the constraints and generate the correct neutrino masses and couplings. For a few benchmark cases we present explicit expressions for the diquark and the leptoquark Yukawa couplings and predict the branching ratios for the rare decays B → φπ and B → φφ. Finally in Sec. V we present our conclusions.

II. COLORED ZEE BABU MODEL
We briefly summarize the main features of the colored Zee Babu model [32,45] that are central to our idea. The model includes a scalar leptoquark S 3L (with lepton number 1) of mass m L and a scalar diquark S D of mass m S transforming as 1 (3, 3, −1/3) and 2 (6, 1, The baryon number of S 3L is taken to be 1/3 whereas S D is assigned 2/3. With this assignment of 1 The choice (3, 1, −1/3) is also possible as it couples neutrinos to down-type quarks but will not explain the R K and R * K anomaly as this scalar couples up-type quarks to charged leptons. 2 Note that if we had chosen the diquark to be (3, 1, −2/3), Y d and, hence, the neutrino mass matrix would be antisymmetric. baryon number, the baryon conservation is automatic and thus the proton decay is forbidden. The lepton number is softly broken through a trilinear term thereby generating Majorana neutrino mass.
With the particle content discussed above, the interaction Lagrangian is given as where α, β = r, b, g are SU (3) c indices, i, j = 1, 2, 3 are generation indices, the diquark coupling matrix, Y ij d , is a symmetric complex matrix whereas the leptoquark coupling matrix, Y ij l , is a general complex matrix. The leptoquark couples to leptons and quarks as √ 2ν iL u jL − √ 2e iL d jL + ν iL d jL + e iL u jl . Note that, in Eq. 3, we can also have additional scalar interaction terms(not relevant to our analysis), such as where Φ is a Higgs doublet. These terms give rise to splitting in the mass of S 3L particles, comprising three states of different electric charges −4/3, −1/3 & 2/3, and thus contribute to the oblique corrections [46]. To avoid that, we assume λ 1,2 = 0 such that all S 3L particles/states have same mass, m L . Along with this, there are quartic and quadratic terms of these scalars. We assume that their coefficients are adjusted such that only the Higgs doublet gets the vev and the potential is bounded from below. neutrino matrix is given as [32,47] where I kl is a loop integral, which in the limit of large leptoquark and diquark masses simplifies to and m d is 3×3 diagonal mass matrix for down-type quarks. Note that we have chosen diagonal bases of the mass matrix for down-type quarks and charged leptons. Hence, to obtain the correct masses of neutrino, we need to diagonalize the mass matrix, M ν by the PMNS matrix U as The standard parametrization is adopted such 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 by the oscillation experiments. However, these phases could be sensitive to the upcoming neutrinoless double beta decay searches.
It should be noted that the mass dimension one parameter, µ, is constrained by demanding the perturbativity of the theory. The trilinear term in the Eq. 3 generates one-loop corrections to leptoquark and diquark masses. These corrections(∆m 2 ) are, in general, proportional to µ 2 16π 2 . Requiring corrections to be smaller than the corresponding masses implies µ ≪ 4πm S/L [47]. As various collider searches, discussed in Sec. V, do not allow the scalar masses to be smaller than 1 TeV, we take µ from 0.1−1 TeV and this choice commensurates with the above constraints.
Having discussed the details of the model, next−we list all the possible constraints, coming from various experiments on leptoquark and diquark coupling matrices.

III. LEPTOQUARKS
• Lepton flavour violation at tree level: Collider searches of leptoquarks indicate that they are heavy. So we can study their low energy effects by writing 4-Fermi operators of two lepton-two quarks. Using Fierz as an effective operator where l and q denote leptons and quarks. These are organized in terms of the four-Fermi effective interactions with normalized dimensionless Wilson coefficients as . In Ref. [48], constraints on such operators have been extensively studied. Keeping in mind that Y ij l should be able to explain a small neutrino mass, following are the most crucial operators related to our work: The µ−e conversion in nuclei sets a bound on the Wilson coefficient of this operator, i.e.
-(µγ µ P L e)(dγ µ P L s): The bound from the decay K • → e + µ − sets a bound on C 1212 The constraint on the K meson decay to pion and neutrinos(ν i ν j ) sets another bound: Apart from this, we have also taken care of all the relevant Wilson coefficients mentioned in Ref. [48].
• Lepton flavour violation radiative decay: The LFV radiative decays l i → l j γ are induced at one loop by the exchange of a leptoquark S 3L with the branching ratio [46] BR where α = e 2 4π , χ µ = 1, and χ τ = 1/5. In the case of a τ lepton, there are two leptonic modes and hadronic modes can be approximated by a single partonic mode(with three colors). Hence there is a factor of 5 difference in µ and the τ -lepton branching ratio. The current experimental bounds [49,50] are • b → sℓ + ℓ − anomalies: As discussed in the Introduction one can perform fits to the b → sℓ + ℓ − data and scenarios in terms of Wilson's coefficients that give a good description of the data. In the above set up, the exchange of the S 3L leptoquark at tree level contributes to the decay b → sℓ + ℓ − , and in particular generates the scenario C µµ 9,NP = −C µµ 10,NP . The effective Hamiltonian describing the decay is parameterized as where O i (µ) are effective operators with Wilson coefficients C i (µ) renormalized at the scale µ. For the model under consideration, only the operators O ℓi 9 = (sγ µ P L b)(l i γ µ ℓ i ) and O ℓi 10 = (sγ µ P L b)(l i γ µ γ 5 ℓ i ) are induced. Using Fierz identity, we obtain the following Wilson coefficients: Assuming new physics only in the muon sector, a model independent analysis on the above operators [17] from the R K , R * K , P ′ 5 and other observables suggests that C µµ 9 (NP) = −0.53 ± 0.08.

IV. DIQUARK
A. Nonleptonic decays and the B → πK puzzle In the Standard Model (SM) the amplitudes for hadronic B decays of the type b → qf f are generated by the following effective Hamiltonian: where the superscript t indicates the internal quark, and f can be a u or c quark. q can be either a d or an s quark depending on whether the decay is a ∆S = 0 or ∆S = −1 process. The operators O q i are defined as where R(L) = 1±γ 5 , and q ′ is summed over u, d, s, candb. O 2 and O 1 are the tree-level and QCD corrected operators, respectively. O 3−6 are the strong gluon induced penguin operators, and operators O 7−10 are due to γ and Z exchange (electroweak penguins) and "box" diagrams at loop level. The Wilson coefficients c f i are defined at the scale µ ≈ m b and have been evaluated to next-to-leading order in QCD. The c t i are the regularization scheme independent values and can be found in Ref. [44].
The diquarks discussed in Sec. II in the context of neutrino mass generation can contribute to the B → πK decays and we can write down the new physics operators that will be generated by a 6 or 3 diquark [51]. In the general case we get the effective Hamiltonian for b quark decays b →d i d j d k as where the superscript d in X d equals 6 or 3 corresponding to the color sextet or the antitriplet diquark. The greek subscripts represent color and the latin subscripts the flavor. We have where the Yukawa Y are symmetric for the sextet diquark and antisymmetric for the antitriplet diquark and we have assumed the same masses for the diquarks.
For b decays of the type b →sss the diquark contribution is tiny as the effective Hamiltonian is proportional to Y d 22 which vanishes for the 3 diquark and is highly suppressed from K and B mixing for the sextet diquark. Similarly the b →ddd transition is proportional to Y d 11 , which is also small. For b → sdd( b →dsd and b →dds) transitions we have the following Hamiltonian: with and We can rewrite the effective Hamiltonian after a color Fierz transformation as The only other unsuppressed transition is b → ssd( b →ssd and b →sds) which has the effective Hamiltonian, with In this case at the meson level we can have the decays B → φπ and the annihilation decays B → φφ. These decays are highly suppressed in the SM and the observance of these decays could signal the presence of diquarks We begin by reviewing the B → πK puzzle. As in Ref. [41] we can analyze the B → πK decays in terms of topological amplitudes. Including only the leading diagrams the B → πK amplitudes become Here, T ′ is the color-allowed tree amplitude, P ′ tc is the gluonic penguin amplitude, and P ′ EW is the color-allowed electroweak penguin amplitude. Furthermore in the SU(3) limit the T ′ and P ′ EW are proportional to each other and so have the same strong phases. Now consider the direct CP asymmetries of B + → π 0 K + and B 0 → π − K + . Such CP asymmetries are generated by the interference of two amplitudes with nonzero relative weak and strong phases.
In both A 0+ and A −+ , T ′ -P ′ tc interference leads to a direct CP asymmetry. On the other hand, in A 0+ , P ′ EW and T ′ have the same strong phase, P ′ EW ∝ T ′ , while P ′ EW and P ′ tc have the same weak phase (= 0), so that P ′ EW does not contribute to the direct CP asymmetry. This means that we expect A CP (B + → π 0 K + ) = A CP (B 0 → π − K + ).
The latest B → πK measurements are shown in Table I. Not only are A CP (B + → π 0 K + ) and A CP (B 0 → π − K + ) not equal, they are of opposite sign! Experimentally, we have (∆A CP ) exp = (12.2 ± 2.2)%. This differs from 0 by 5.5σ. This is the naive B → πK puzzle.  for the four B → πK decay modes. The data are taken from Ref. [52].

C. Model-independent new physics formalism
In the general approach of Refs. [53,54], the NP operators that contribute to the B → πK amplitudes take the form O ij,q NP ∼sΓ i bqΓ j q (q = u, d), where Γ i,j represents Lorentz structures, and color indices are suppressed. The NP contributions to B → πK are encoded in the matrix elements πK| O ij,q NP |B . In general, each matrix element has its own NP weak and strong phases.
Note that the strong phases are basically generated by QCD rescattering from diagrams with the same CKM matrix elements. One can argue that the strong phase of T ′ is expected to be very small since it is due to selfrescattering. For the same reason, all NP strong phases are also small, and can be neglected. In this case, many NP matrix elements can be combined into a single NP amplitude, with a single weak phase: Here the strong phase is zero. There are two classes of such NP amplitudes, differing only in their color structure: They are denoted A ′,q e iΦ ′ q and A ′C,q e iΦ ′C q , respectively [54]. Here, Φ ′ q and Φ ′C q are the NP weak phases. In general, A ′,q = A ′C,q and Φ ′ q = Φ ′C q . Note that, despite the "color-suppressed" index C, the matrix elements A ′C,q e iΦ ′C q are not necessarily smaller than A ′,q e iΦ ′ q .
There are therefore four NP matrix elements that contribute to B → πK decays. However, only three combinations appear in the amplitudes: [54]. The B → πK amplitudes can now be written in terms of the SM diagrams and these NP matrix elements. Here we neglect the small SM diagram P ′ uc but include the color-suppressed amplitudes: In Ref. [55], a different set of NP operators is defined: In this case we have  f. and best-fit values of unknown parameters for the Diquark model where the fit 1 has X 6 = X 3 , and fit 2 has X 3 = 0. Constraints: B → πK data, measurements of β and γ, |C ′ /T ′ | = 0.2, |P ′C EW,N P /P ′ EW,N P | = 0.3 (fit 1), and |P ′C EW,N P /P ′ EW,N P | = 1 (fit 2).
We consider two models, the first with This leads to P ′C EW,N P /P ′ EW,N P = 1 3 with both amplitudes having the same weak phase.
The second model has This leads to P ′C EW,N P /P ′ EW,N P = 1, again with both amplitudes having the same weak phase. A χ 2 fit for the new physics within this scenario is performed to determine the parameters of the model. The procedure for determining such a fit is as follows. We define the function where O exp and ∆O exp are the experimentally determined quantities with their associated uncertainties, respectively, as listed in Table I. O th are determined from the model and are thus functions of the unknown parameters. The goal from here is to find the values of the parameters that minimize χ 2 . There are many programs available to accomplish this, one of the most widely used is MINUIT [56], which is used here. The goodness of the fit is determined by the value of χ 2 at the minimum and the number of degrees of freedom in the fit. The degrees of freedom are the number of constraints included in the fit minus the number of parameters that are fitted. In this case the number of constraints is 13: the B → πK data, the independent measurements of β and γ, and the constraints on |C ′ /T ′ | and |P ′C EW,N P /P ′ EW,N P |. The number of parameters is nine and we have that the number of degrees of freedom are four. A "good" fit is one where χ 2 min ≈ d.o.f., but a better measure is the p value which gives the probability that the model tested adequately describes the observations.
The results of the fit for this case are shown in Table II. Here the p value is 44% for X 6 = X 3 , and 43% for X 3 = 0, which is not bad (and is far better than that of the SM).
The SM T ′ diagram involves the tree-level decayb →ūW + * (→ us = K + ). The NP P ′ EW,N P diagram looks very similar and is expressed relative to the T ′ diagram. Within factorization, the SM and NP diagrams involve (0)f π , respectively, where F B→K,π 0 (0) are form factors and f π,K are decay constants. The hadronic factors are similar in size: |A Kπ /A πK | = 0.9 ± 0.1 [44]. Taking central values for X 6 = X 3 , we have [41] For X 3 = 0 we obtain Both models give similar fits and in Fig. 2 we show the allowed regions of the diquark couplings within a 1σ range for the first model.

D. Neutral meson Mixing
Diquarks, in spite of being charged, through their coupling to the same generation quarks can mediate the mixing between neutral mesons at tree level. Following the convention in [57], the mixing can be depicted as the six dimension operator: The 90 % C.L. bounds on the corresponding Wilson coefficients [57] is then given as

V. NUMERICAL ANALYSIS AND DISCUSSION
Before we present the results, we discuss the bounds on the scalar masses obtained from collider experiments.
The collider experiments provide direct limits on the leptoquark mass when they decay to leptons and quarks in the final state. There are many studies in the literature where different signatures have been discussed [32,58,59].
The leptoquarks can be pair produced from gg and qq as initial state or singly produced at hadron colliders via g + q → S 3L + lepton. Recent studies at ATLAS [60] and CMS [61] with 13 TeV data puts a bound on the scalar leptoquark mass, m L > 1, 1.2 (ATLAS), 0.9 (CMS) TeV when decay to u e, c µ, and t τ with 100% branching fraction, respectively, at 95 % C.L. The previous results [62,63]  As for the leptoquark case, Y 2i l couplings(real) are generated randomly in the range [10 −5 : 1]. With the obtained sets of couplings, we calculate the strength of remaining leptoquark couplings, for randomly generated LQ mass, from Eq. 4 to get the correct neutrino masses. The symmetric neutrino mass matrix in the Eq. 4 represents six independent equations as six independent parameters (given in Table III) that are obtained from the neutrino oscillation experiments. Throughout the analysis, we have kept Majorana phases to be 0, and have employed the 2σ ranges for the neutrino mixing parameters for normal hierarchy from Refs. [65,66]. Finally, those sets of LQ couplings  We compare our results for leptoquark coupling with the results given in [33] and [67] and find them consistent. A few benchmark points (BP) are given in Appendix A. For these BP, we present branching ratios for the rare decays in Table IV following the calculations in Ref. [51]. The branching ratios are rather small and it will be difficult to observe these decays in ongoing experiments. Our analysis shows that the B anomalies and the neutrino masses can all be accommodated in a consistent framework.

VI. CONCLUSION
In conclusion we have discussed a unified framework to provide solutions to three problems. They are the anomalies in b → sµ + µ − measurements, nonleptonic B → πK decays, and the issue of generating neutrino masses and mixing.
Our framework contained a scalar triplet leptoquark, a scalar color sextext diquark, and also, possibly, a color antitriplet diquark. We considered several low energy as well as collider bounds on the leptoquark, diquark couplings and masses. For the leptoquarks these low energy observables included the b → sℓ + ℓ − measurements. The solutions to the B → πK puzzle provided constraints on products of the diquark Yukawa couplings. We then checked that the correct neutrino masses and mixings were reproduced with the allowed couplings of the leptoquarks and diquarks.
We also predicted the branching ratios for a few rare B decays whose observations could signal the existence of diquarks. However, we found the branching ratios of these decays to be unobservably small. Because S αβ is symmetric/antisymmetric there is an additional factor of 2. In other words S 12 can contract with S 12 and S 21 .

Appendix B: Benchmark Points
Here we give the benchmark points satisfying the B anomalies observations and explaining the neutrino mass.