Combined explanations of $(g-2)_{\mu}$, $R_{D^{(\ast)}}$, $R_{K^{(\ast)}}$ anomalies in a two-loop radiative neutrino mass model

Motivated by the long-standing tension in the muon anomalous magnetic moment (AMM) and persistent observations of B-physics anomalies in $R_{D^{(\ast)}}$ and $R_{K^{(\ast)}}$ ratios, we construct a simple two-loop radiative neutrino mass model, and propose a combined explanations of all these apparently disjoint phenomena within this framework. Our proposed model consists of two scalar leptoquarks (LQs), a $SU(2)_L$ singlet $S_1\sim (\overline{3},1,1/3)$ and a $SU(2)_L$ triplet $S_3\sim (\overline{3},3,1/3)$ to accommodate $R_{D^{(\ast)}}$ and $R_{K^{(\ast)}}$ anomalies, respectively. The muon receives chirality-enhanced contribution towards its $g-2$ due to the presence of $S_1$ LQ that accounts for the observed deviation from the Standard Model prediction. Furthermore, we introduce a $SU(2)_L$ singlet scalar diquark $\omega\sim (\overline{6},1,2/3)$, which is necessary to break lepton number and generate neutrino mass radiatively with the aid of $S_1$ and $S_3$ LQs. We perform a detailed phenomenological analysis of this set-up and demonstrate its viability by providing benchmark points where a fit to the neutrino oscillation data together with proper explanations of the muon AMM puzzle and flavor anomalies are accomplished while simultaneously meeting all other flavor violation and collider bounds.


Introduction
In the Standard Model (SM), contributions to the anomalous magnetic moment (AMM) of the muon a µ , arising from loop corrections [1] are calculated with excellent accuracy. On top of that since experiments determine this quantity to high precision, any deviation from the theory prediction directly points towards physics beyond the SM (BSM). In fact, there is a long-standing discrepancy between the theoretical computations [2][3][4] and its measured value [5], corresponding to a 3.7σ anomaly. In the coming days, the Muon g − 2 Collaboration [6] at Fermilab is expected to announce their result, which further motivates our investigation of the possible NP explanation of this anomaly.
Over the last two decades, various mechanisms are proposed to account for this deviation.
Among them the effects of scalar leptoquarks (LQs) on a µ are studied extensively, see for example Refs. [7][8][9][10] for single LQ solution to (g − 2) µ . LQ extensions of the SM has gained a lot of attention recently, due to their ability in accommodating the persistent tensions observed in the lepton flavor universality violating B meson decays, particularly in the R K ( * ) and R D ( * ) ratios (see for example Refs.  for both scalar and vector LQs explanations). These anomalies include flavor changing neutral current b → s, as well as flavor changing charged current b → c transitions, which we briefly summarize below.
Recent measurements have observed notable digressions from the SM predictions in the following two ratios associated with neutral current transition: Theory predictions of these ratios are: R SM K = 1.0003 ± 0.0001 [80], R SM K * = 1.00 ± 0.01 [81].
Even though these values are in harmony with the SM, results from the LHCb show significant deviations compared to theory predictions, 110 −0.070 ± 0.024, 0.045 GeV 2 < q 2 < 1.1 GeV 2 [84], 0.685 +0. 113 −0.069 ± 0.047, 1.1 GeV 2 < q 2 < 6.0 GeV 2 [84]. (6) These measured values in both low and high q 2 bins point towards 2.5σ deviation from SM values. Discrepancies observed in the R K and R K * ratios have gained much curiosity in the theory community due to their trustable theory predictions, since hadronic uncertainties cancel out in these ratios.
The outstanding tension of the muon AMM together with the large deviations measured in the lepton flavor universality violating decays of the B mesons clearly indicate the existence of new physics beyond the SM. As already aforementioned, scalar leptoquarks are the prime candidates in resolving these observed anomalies. However, a single scalar LQ cannot accommodate for three of these anomalies simultaneously. First, we identify the pair of LQs that can do our desired job. For a TeV scale LQ, a large enough contribution is required to account for ∆a µ data, which can be provided if both the left-handed and right-handed chiral couplings of the LQ are present [10]. This requirement is satisfied by only two scalar LQs, S 1 ∼ (3, 1, 1/3) and R 2 ∼ (3, 2, 7/6). It is interesting to realize that either of these two LQs can accommodate anomalies in the R D and R D * ratios at the tree-level (see for example Ref. [56]). On the other hand, S 3 ∼ (3, 3, 1/3) is the only scalar LQ that can correctly incorporate R K and R K * anomalies at the tree-level (see for example Ref. [56]).
By following the above discussion, in this work, we postulate that the NP beyond the SM contains S 1 and S 3 LQs. With these in hand, one must ask the obvious question: how to give mass to the neutrinos 1 ? The reason for this is, even though neutrinos remain massless in the SM, observations of neutrino oscillations are securely established by a number of experiments [101][102][103][104][105][106][107]. Hence, any BSM construction is obliged to explain the origin of neutrinos masses and mixings. It gives rise to a more appealing scenario if the BSM states introduced to resolve these tensions also participate in neutrino mass generation mechanism 2 . Since proper explanations of the above-mentioned anomalies demand TeV LQs, it is evident that the only natural choice to generate neutrino mass is via quantum corrections [108][109][110][111][112][113][114]. However, with just S 1 and S 3 LQs added to the SM, neutrinos cannot get mass 3 . We must introduce one more BSM particle in the theory. One obvious and simple choice is to extend the scalar sector by a color sextet diquark 4 ω ∼ (6, 1, 2/3), which is a singlet under the SU (2) L . Addition of this scalar diquark (DQ) breaks the lepton number by two units, and Majonara mass for the neutrinos are then generated at the two-loop level, in which all three BSM particles run through the loop.
In a nutshell, we propose a framework in which the neutrino mass, the muon anomalous magnetic moment puzzle, and B-physics anomalies in the R D ( * ) , R K ( * ) ratios have a common origin. We perform a comprehensive phenomenological analysis of this set-up and discuss the feasibility of interpreting these anomalies as well as explaining the neutrino oscillation data.
In the next section (Sec. 2), we introduce the model and then discuss how to ameliorate these anomalies in Sec. 3. Relevant experimental constraints on the model parameters are detailed in Sec. 4. We present the results in Sec. 5 and finally conclude in Sec. 6.
1 Instead of S 1 , if R 2 is used in association with S 3 , neutrino mass generation and reconciling B-physics anomalies are discussed in Refs. [70,78]. 2 See for example Refs. [18,21,33,38,53,69,70,78,79] that unify neutrino mass generation mechanism with B-physics anomalies. 3 Extension of the SM with S 1 and S 3 LQs was considered in Ref. [75] without addressing the question of neutrino mass generation. On the other hand, in Ref. [71], vectorlike-quarks ∼ (3, 2, −5/6) was introduced in addition to S 1 and S 3 LQs to give neutrinos non-zero masses. 4 Ref. [115] first proposed neutrino mass generation at the two-loop by introducing S 1 and ω. This model is then analysed in more details and collider implications of these new colored states are studied in Ref. [116].
Ref. [38] considered the scenario of utilizing S 3 instead of S 1 in neutrino mass generation and to incorporate only R K ( * ) anomaly. Furthermore, Ref. [69] had the same particle content as that of Ref. [38] and their work focused on explaining R K ( * ) and B → Kπ anomalies. None of these frameworks can simultaneously explain the tensions in the R D ( * ) and a µ , which we attempt to achieve in this work.
The existence of S 1 LQ can account for the anomaly observed in the muon AMM a µ .
Furthermore, both these LQs accompanied by the DQ ω participate in generating neutrino mass radiatively at the two-loop level, as shown in Fig. 1. As already aforementioned, existence of the DQ is required to break lepton number by two units, and provide mass to the neutrinos. Hence, in our model neutrinos are Majorana like fermions. The Yakawa couplings associated to the LQs are given as follows [117]: as usual, here Q and L are the left-handed quark and lepton doublets of SU (2) L , and d R , u R , and R are the right-handed down-type quark, up-type quark, and charged lepton, respectively, which are all singlets of SU (2) L . Here σ a (a = 1, 2, 3) are the Pauli matrices, and {i, j} are flavor indices. Moreover, S a 3 are the components of S 3 in the SU (2) L space. In the above Lagrangian, we have omitted the S 1,3 couplings to diquarks to ensure proton stability. The Yukawa couplings y L , y R and y S are a priori arbitrary 3 × 3 matrices in the flavor space.
To calculate the flavor observables, it is convenient to write the above Lagrangian in the charged fermion mass eigenbasis, for which we make the following transformations of the fermion fields: Here U and V are the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) and Cabibbo-Kobayashi-Maskawa (CKM) matrices, respectively. Moreover, following the notation of Ref. [33] we have definedν that represents the neutrino weak-eigenstate. With these, the Lagrangian takes the following form: In the above set of Lagrangian terms, the Yukawa couplings of the DQ scalars are also summarized, which will be required for neutrino mass generation. Note that y ω is a 3 × 3 symmetric matrix in the flavor space.
Before proceeding any further, here we clarify few assumptions that we make. First, we assume all Yukawa couplings to be real for simplicity. The electroweak (EW) symmetry breaking will split the masses of the three component fields that belong to the triplet LQ.
However, splittings among different components are highly constrained by EW precision measurements, this is why we chose them to be degenerate in mass. Subsequently, we ignore any mixing between the S and ω, respectively. 3 . These two-loop neutrino mass diagrams utilize the following cubic coupling terms in the scalar potential: With these, the neutrino mass formula has the following form [116]: Here m d = diag{m d , m s , m b } is the diagonal down-quark mass matrix. In this formula there are two terms, one for p = 1 for which we have µ p = µ 1 , y p = y L , and the second for p = 3 that corresponds to µ p = µ 3 , y p = y S . Since the down-quark masses are very small compared to the LQ and DQ masses, the loop integrals have the following simple and generation independent form [118]: Here M p = M 1 (M 3 ) for p = 1 (p = 3). Since the mass generation occurs at the twoloop level, TeV scale BSM states running in the loop naturally provide tiny masses to the neutrinos without requiring the Yukawa couplings to be abnormally small. In fact as we will show, Yukawa couplings of order 0.01 − 1 are the prerequisite for concurrent explanations of B-physics anomalies, muon AMM, as well as neutrino oscillation data.
3 Resolving Anomalies When considering LQ solutions to lepton AMMs, it is well known that relevant contributions can only be provided by a non-chiral LQs [10] as discussed above. Among the only two nonchiral LQs R 2 and S 1 , the latter is present in our set-up. Even though both S 1 and S 3 can in principle contribute to a within our scenario, only effects coming from S 1 are important due to chiral enhancement. The dominant one-loop contributions to charged lepton AMM are presented in Fig. 2. The effective Lagrangian from which a is calculated can be written as [10,119,120]: where the field strength tensor is defined as Here the contributions σ L,R can be computed from the effective Lagrangian that leads to → γ, which is given below [10,119,120], To a very good approximation we find the corresponding contributions relevant to our study have the following expressions (by setting = ): where we have assumed all couplings to be real. This leads to the following expression for the muon AMM arising dominantly from S 1 LQ: Here we have used the color factor N c = 3, V tb = 1, and x t = m 2 t /M 2 1 . Due to the topquark mass insertion inside the loops as shown in Fig. 2, the observed enhanced value of the muon magnetic moment can be naturally incorporated within this framework for a TeV scale leptoquark.

R K and R K *
It is remarkable that S 3 is the only scalar LQ that can simultaneously account for R K < R SM K and R K * < R SM K * at tree-level. Processes of the form B → K ( * ) + − can be described by the following effective Hamiltonian where the effective operators are given by After integrating out the heavy leptoquark and combining the Yukawa part of the Lagrangian associated to S 3 as given in Eq. (13) with the above effective Hamiltonian lead to the following purely vector Wilson coefficients at the LQ mass scale: By assuming the NP coupling to electrons is negligible (leading to = = µ), the observed values of the R K and R K * ratios then can be explained with C 22 9,10 < 0. Wilson coefficients of this type are generated within our framework by the S 3 LQ couplings to muons over the electrons as depicted in Fig Discrepancies are also founds in few other observables related to neutral current processes. For example the most significant departure has been found in the angular observable P 5 in the B → K * µµ decay [122][123][124]. Unlike R K ( * ) , these additional quantities suffer from hadronic uncertainties and we do not include them in our analysis.

R D and R D *
As for the charged current process b → cτ ν that is responsible for B meson decays B → Dτ ν and B → D * τ ν get contributions from both S 1 and S 3 LQs at the tree-level. Feynman diagrams that lead to such processes are shown in Fig. 3 (right diagram). Processes of these types can be described by considering the following effective Hamiltonian: where in the SM C SM V = 1. In the above effective Hamiltonian, both S 1 and S 3 contribute to the vector Wilson coefficient, whereas only S 1 participates in the scalar and tensor Wilson coefficients, which at the LQ mass scale have the following forms: We will focus on scenarios where dominant coefficients are the ones with i = 3, which corresponds to lepton flavor conservation [12,17,33,[125][126][127][128]. Then the expressions of the R D and R D * ratios are given by [129]: In these formulas, the Wilson coefficients are given at the low scale µ = m b .
In addition to R D ( * ) , there are a number of observables associated to the charged current processes that indicate disagreements to some extent when compared to the SM values.
Such as the ratio R J/ψ of the tauonic mode to the muonic mode for B → J/ψ ν [130], the longitudinal polarization of the D * denoted by f D * L [131], and polarization asymmetry in the longitudinal direction of the tau in the D * mode denoted by P * τ [94]. These observables have comparatively large error bars, hence we focus only on explaining R D ( * ) .

Synopsis
Here we discuss the textures of the Yukawa coupling matrices required for a combined explanations of the aforementioned phenomena that we want to achieve in this work. First note that three Yukawa coupling matrices, y L,S,ω enter in the neutrino mass formula given in Eq. (16). Among them, y L participates in explaining both the muon AMM and R D ( * ) anomalies, whereas y S takes part in incorporating R K ( * ) and R D ( * ) ratios. The only Yukawa coupling matrix, y R that does not contribute to neutrino masses, however plays significant role in resolving tensions in a µ and R D ( * ) .
As for the neutrinos, two mass squared differences and three mixing angles have been measured with great accuracy. Even though the hierarchical pattern, whether normal or- is not yet known, inverted ordering is less favored by the data. Hence in this work, we stick to normal ordering for the neutrinos. In addition to masses and mixings, if neutrinos are Majorana fermions, then there are three more physical quantities exists in the neutrino sector. One of them is the Dirac phase, and rest two are Majorana phases. Dirac phase, which has not been measured yet directly, currently has large uncertainty associated to it [132], furthermore we have no clue about the range of the Majorana phases. In this work, we take all parameters to be real, and do not focus on predicting these phases. An overview of the most recent global fit [133] to the neutrino oscillation data are given as follows: Following the aforementioned discussions on reconciling these anomalies along with neutrino oscillation data, we adopt the following form of the Yukawa coupling matrices: The entries in blue plays role in ∆a µ , the entries in green enters into the R K ( * ) expressions, and the couplings in red contribute to R D ( * ) ratios. Additionally the entries in black are introduced to get consistent fit to the neutrino masses and mixing angles. A few comments are in order regarding the choice of the above Yukawa couplings matrices. Since explanations to the B meson decay anomalies demand the existence of most of the entries in the lower 2 × 2 blocks, it is a natural choice to populate y ω matrix in the same lower 2 × 2 block.
We intentionally do not introduce any couplings with the first generation of quarks in the above matrices, since these couplings are severely constrained by many different experiments.
However, it is easily understood that filling out all entries in the lower 2 × 2 block are not sufficient to give a realistic fit to the neutrino data. Introducing y ω 21 or y ω 31 term does not change the above conclusion either. This leaves us with four different minimal options, considering one non-vanishing term from the set {y S 21 , y S 31 , y L 21 , y L 31 }. Instead of exploring all such possibilities, we fix y S 31 = 0 for the rest of the analysis. In the next section, we elucidate the experimental constraints on the aforementioned non-zero Yukawa couplings.
Before closing this section, here we briefly discuss the loop corrections and running of the Wilson coefficients. The QCD corrections to the matching on 2-quark-2-lepton operators mediating semileptonic B decays have been recently computed in Ref. [134]. This correction leads to a shift of the Wilson coefficients Eqs. (25), (27), and (28) that are, These QCD corrections enhance the contributions to about 10% [134] which definitely favor towards the explanations of B meson decay anomalies. Furthermore, to evaluate the abovementioned flavor ratios, we run these operators to the bottom-quark mass scale at which the relevant form factors are calculated. We use the Flavio package [135] to do this running (see also Ref. [136]) and find the following relations between the two different scales: In this calculation, we have fixed M LQ = 1200 GeV, and for the bottom-quark mass m b = 4.18 GeV is used. The relation between the scalar and the tensor Wilson coefficients also gets modified at the low scale, which we find to be C S (µ = m b ) = −7.63 C T (µ = m b ).

Correlated Observables
In the previous sections, we have discussed the NP contributions to the muon AMM and R K * , R D * flavor ratios, and the neutrino mass generation mechanism is introduced in Sec.
2. Accommodating these significant deviations from the theory predictions lead the way to various flavor violating processes that are severely constrained by experimental data. In this section, we consider all such relevant processes and discuss the associated constraints on the model parameters.

→ γ Processes
The effective Lagrangian leading to radiative decays of the charged leptons → γ is given in Eq. (20). Although both LQs mediate these dangerous processes, S 1 mediated τ → µγ receives chirality-enhanced effect from top-quark, which is the strongest constraint within our model. The branching ratios associated to these process are calculated by the following formula [119]: where τ is the lifetime of the initial state lepton and we derive the following expressions of these σ L,R originating from S 1 and S 3 LQs [119,120]: Here as before x q = m 2 q /M 2 LQ . However, for q = t, the replacement of m q → µ LQ inside the log-function needs to be made for consistency, see Ref. [75] for details. In the above formulas terms proportional to y R (y R ) * are not shown, since they vanish for our choice of the Yukawa textures. In the following we summarize the current experimental limits on these processes [137,138]: Br (τ → µγ) < 4.4 × 10 −8 .

→ Processes
The interaction terms that lead to → γ also generate the rare lepton flavor violating decays → . LQs present in our set-up induce these processes at the one-loop level. Decays of these types proceed via penguin-diagrams with Z and γ exchanges, and via box-diagrams with quarks and LQs inside the loops. The corresponding box-diagram contributions are always somewhat smaller than the penguin-diagrams, hence we omit those terms. Branching ratios of such decay channels can be written as [139][140][141][142][143]: A slight modification of the above expression is required when there are two different lepton flavors in the final state [142,143]: Photon (Z-boson) penguin-diagrams are encoded in the T 1L,1R and T 2L,2R (Z L,R ) terms. We derive the following expressions of these terms from S 1 and S 3 LQs: Here we have defined: , and g d S 3 = 4 sin 2 θ W /3. Moreover, θ W is the Weinberg angle. The current experimental bounds of these processes are quoted below [144,145]: Br(τ ± → µ ± e + e − ) < 1.5 × 10 −8 .

Z → Processes
The Z-boson decays to leptons receive contributions from the LQs that constraint the Yukawa couplings. These processes are explained with the following effective Lagrangian: Here g is the SU (2) L gauge coupling. These dimensionless couplings g ij are very accurately measured at the LEP [146] that provide stringent constraints on the associated Yukawa couplings for a fixed LQ mass. NP contributions to these dimensionless couplings can be expressed as follows [147]: When calculating δg L we have defined Similarly while calculating δg R we make the replacements w u ij = y R ij , The results from the LEP collaboration [148] provide the following limits on the NP contributions: Furthermore, the branching ratio for the processes Z → are given by [147]: . Both LEP and LHC results put limits on these branching ratios which are [149][150][151] summarized below: As for the neutrinos, Z-decays of the form Z → νν also receive contributions form LQs that are parametrized by, Above we have also collected the accurately measured experimental value of this effective number of neutrinos.

µ − e Conversion
With the choice of the Yukawa coupling matrices given in Eq. (33), S 3 LQ mediates µ − e conversion in nuclei at the tree-level in our model. This conversion rate can be calculated from the following formula [78,120,152]: Γ capture (Z) is the total capture rate for a nucleus with atomic number Z, which is 13. µ ). The current sensitivity implies CR(µ − e) < 7 × 10 −13 [153], whereas the future projected sensitivity is expected to make almost four orders of magnitude improvement over the current limit CR(µ → e) < 10 −16 [154][155][156][157][158][159][160].

P 0 → − +
For the explanations of the R K ( * ) ratios along with neutrino oscillation data, the NP contributions to the O 9,10 operators need to be large. The associated Wilson coefficients C 9,10 as given in Eq. (25) then lead to interesting pseudoscalar meson decays via b → sµ + µ − , b → sµ + τ − , and b → sτ + τ − . The decay width of the process P 0 → − + can be written as [161]: For our scenario B s is the only relevant meson, which corresponds to q = t, j = b, i = s in the above formula, and the function η is defined as: The experimental limits on these processes are given below [162][163][164]: Among these, only B s → µ + µ − decay mode has been observed, which is in good agreement with the SM prediction [165], Br(B s → µ + µ − ) SM = (3.65 ± 0.23) × 10 −9 .
Associated to b → sµ + τ − transition there is another important constraint that comes from B → K decay that has the following branching ratio [166]: with the following experimental bound on this process [167]:

B → K ( * ) νν
Both S 1 and S 3 LQs can induce B → K ( * ) νν decay at the tree-level via d k → d j νν processes.
The Wilson coefficients responsible for such decays associated to b → s transition are: Then following [15], the branching ratio for B → K ( * ) νν can be expressed as: this ratio is normalized to SM, where C SM L = −1.47/ sin 2 θ W . The Belle collaboration [168] limits these ratios to be R νν K < 3.9 and R νν K * < 2.7.

B c → τ ν
The same Wilson coefficients that explain R D ( * ) in our framework also lead to B c → τ ν decay. The associated branching ratio that depends on the vector and the scalar Wilson coefficients can be written as [129,169]: The lifetime of B c has not been measured in the experiments yet. Hence, this quantity needs to be compared with the theoretical calculations [156,[170][171][172][173]. By carrying out such calculations in Ref. [174] and Ref. [175], their results advocate that the NP contributions to this decay must be Br(B c → τ ν) ≤ 10% and Br(B c → τ ν) ≤ 30%, respectively. On the other hand, as argued in Refs. [129,176], these calculations suffer from theoretical uncertainties, and suggested a conservative limit of Br(B c → τ ν) ≤ 60%. It is interesting to note that R 2 LQ explanations to R D ( * ) demands much larger values [78] of this branching ratio, hence such explanations can in principle be ruled out by reducing the corresponding uncertainties in future. On the contrary, the observation of R D ( * ) can be properly accommodated via S 1 LQ with smaller values of this branching ratio [75].

τ → P 0 Decays
In the SM tau lepton decays into mesons and lighter leptons are not allowed. However, these lepton flavor violating decays can be significant in the presence of leptoquarks that provide strong constraints on the Yukawa couplings. Tau lepton decay width for τ → P 0 process can be written as follows [143]: where λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2ac − 2bc. Within our scenario, the related processes we need to take into account are for P = φ, η, η . As for the meson form factors we take the number quoted in Ref. [26], and their masses are taken from Ref. [146]. From hereafter, we will neglect the mass of the lighter charged lepton. With this assumption, the only relevant terms that enter in Eq. (85) are: Current bounds on these processes are [146], Concerning the LQs, both S 1 and S 3 contribute to meson-antimeson mixing. This NP contribution to B 0 s − B 0 s mixing can be described by the effective Lagrangian given below [177]: Here the SM part is C SM [178] and the NP contribution at the heavy scale (Λ) is given by [75,120,177,179]: Here we neglect the evolution of C N P 1 from high scale to the m w scale. Then the mass difference is given by: The SM prediction is ∆m SM Bs = (18.3 ± 2.7) × 10 12 s −1 [180,181]. This mass difference has been measured in the experiments [146,182] with great accuracy, which is given by:

Numerical Analysis and Discussion
In this section we perform a numerical analysis of the proposed model to demonstrate how to reconcile neutrino oscillation data with anomalies in the B meson decays and the muon anomalous magnetic moment.
First we recall that for simplicity, we consider all parameters of this theory to be real.
Extension to the general case with complex Yukawa couplings is straightforward. In this CPconserving scenario, we adopt the Wolfenstein parametrization [183] for the CKM matrix: and take values of the mixing parameters λ=0.2248, A=0.8235, ρ=0.1569 [146].  This corresponds to m 0 = 1.95 × 10 −10 µ, then µ can be fixed from one of the two measured neutrino mass squared differences.
As for the neutrinos alone it is trivial to get a fit to the data from the above mass matrix formula. However, as elaborated in the previous sections, the explanations of the muon The type of Yukawa coupling textures that we consider in this work is already introduced in Eq. (33). As we have discussed in Sec. 3, in a scenario with only entries in the lower 2 × 2 blocks for all the matrices does not lead to realistic neutrino fit. It is trivial to understand that two of the three mixings angles θ 12 and θ 13 would remain zero in this case. As we have argued, to alleviate this issue, one needs to consider at least one non-vanishing term among {y S 21 , y S 31 , y S 21 , y S 31 }, and we have made an ad hoc choice of y S 31 = 0 just for demonstration. An immediate consequence is that non-zero y S 31 leads to µ → e conversion in the nuclei. Hence, neutrino oscillations are directly linked to lepton flavor violating processes. Correlations among these quantities are depicted in Fig. 4 by randomly varying the relevant Yukawa couplings.
Since within our scenario, both the vector and scalar-tensor Wilson coefficients take part in explaining the R D ( * ) ratios, significant new physics contributions to B 0 s − B 0 s mixing, as Similarly, blue (yellow) dots the case where NP contribution to B 0 s − B 0 s mixing is in between 10% and 20% (20% and 50%). From this plot, it is clear that a fit to both R D and R D * to their experimental central values require more than 10% contribution to B 0 s − B 0 s mixing. In the same figure, the plot on the right shows the interrelationship between R D and the branching ratio for B c → τ ν. As can be seen from this plot that correct values of R D and R D * ratios can be reproduced within this set-up even with Br(B c → τ ν) < 10%, which is unlike the scenarios when R 2 LQ is employed to explain B meson decay anomalies in the charged current processes that demands large branching ratio of this process (see for example Ref. [78]). Another immediate difference between utilizing R 2 and S 1 that we point out here is, even though in our scenario for certain choices of parameters, NP contributions to Z → τ L(R) τ L(R) can be large, consistent fits can be obtained where these relevant contributions are small (see Table I). However, when S 1 is replaced with R 2 LQ, NP effects on Z → τ L(R) τ L(R) decays are usually significant that puts strong restrictions on the upper limit on the associated Yukawa couplings (see for example Ref. [78]).
The essential parameters that describe the muon AMM and B-physics anomalies, as well as neutrino oscillation data unavoidably lead to charged lepton and meson decays. In our set-up, the tau decays to a muon and a photon is the most constraining process. In fact, as long as τ → µγ decay limit is satisfied, τ → µµµ processes are under control.
Interconnections among tau decays to lighter leptons and a photon, as well as its decays to meson and lighter leptons are presented in Fig. 6 by varying the relevant Yukawa couplings. Further correlations among different meson decay modes are portrayed in Fig. 7.
From our detailed numerical analysis we find that points that satisfy all fit requirements, branching fractions for τ → µγ, τ → φµ and τ → ηµ are always very close to the current experiment upper limits (see Table I). As for the lepton flavor violating B meson decays, B s → τ µ is expected to be within one or two orders below the current experimental limit, whereas for decays of the form B → Kτ µ, the expected branching ratios are just one order smaller than the current bounds (see Table I). Furthermore, NP contributions to the branching ratios of B s → τ τ and B s → µµ are about hundred times enhanced and suppressed, respectively compared to the SM predictions (see Table I). Some of these enhanced effects, such as in B → Kτ µ and τ → φµ can be tested soon by LHCb and Belle-II collaborations.
For illustration purpose, we also provide concrete benchmark points that incorporate neutrino masses and mixings, as well as accommodate anomalies in the a µ , R K ( * ) , R D ( * ) , and simultaneously satisfy all experimental constraints. Two such benchmark points are given below:  Table I.

LHC Bounds
In this sub-section, we briefly discuss the collider bounds. There exists dedicated direct search for LQs at the LHC that provide strong bounds on the masses of the LQs. From our numerical inspection the typical types of solutions that we get are of similar forms as the benchmark points presented in Eqs. (99) and (100). This is why it is sufficient to discuss the representative bounds associated with these benchmark points. Since the Yukawa couplings are not too large in our scenario, hence the main LHC bounds are coming from the QCD driven LQ pair-production. Neglecting the t-channel contributions, this corresponds to two different production mechanism: gluon-gluon fusion (gg → LQLQ), and quark-antiquark annihilation (qq → LQLQ). Once produced, each LQ will decay into a quark and a lepton, and the bounds on these LQ masses highly depend on the branching fractions to different decay modes.
For illustration let us take for example BM-II of Eq. (100) to derive these bounds. The decay modes of the LQs are then given by: here numbers inside the parentheses are the associated Yukawa couplings responsible for these decay modes. The bounds on the LQ masses for these decay modes from LHC searcher are given as follows: LQLQ → jjνν : 980 GeV (635 GeV); @35.9 f b −1 [184], .
Here the current limits on LQ masses are shown for 100% (50%) branching ratios. "j" represents a jet that could be any light quark, for example u, d, s, c. Moreover, LHC luminosity for each search along with the experimental references are shown for each decay modes.
Even though LQLQ → ttµµ decay mode has the largest bound on LQ mass of 1420 GeV, within our scenario, the corresponding branching ratio is negligibly small, leading to a much lower mass bound. Consequently, the chosen leptoquark mass of M LQ = 1200 GeV for our numerical analysis safely satisfies all collider bounds.
As for the diquark, LHC searches for dijets in the final state. Diquark mass smaller than 6 TeV is ruled out by recent collider studies [189]. It should be pointed out that this limit was derived for diquarks that has couplings to up-quarks, which in our case only couples to down-quarks. Consequently, the lower bound on the mass is expected to be somewhat smaller. Not to mention, the limits largely depend on the associated branching ratios. In this work, for simplicity, we assume its mass to be much heavier compared to the LQs such that all collider, as well as other experimental constraints are automatically satisfied. For example, among all processes mediated by the DQ, the most dangerous constraint comes from its contribution to B 0 s − B 0 s mixing. Following Ref. [190] we find very strong bounds on the Yukawa coupling that are given by: This for M ω = 10 TeV leads to |(y ω 22 ) * y ω 33 | < 0.2. As can be easily verified from the benchmark points provided in Eqs. (99) and (100), the types of solutions we achieve meet all requirements.

Conclusion
In this work we have explored the possibility that the neutrino mass, the long-standing tension in the muon anomalous magnetic moment, and persistent observations of B-physics anomalies in the R D ( * ) , R K ( * ) ratios have a common origin. Our proposal is a simple extension of the Standard Model that consists of two scalar leptoquarks S 1 ∼ (3, 1, 1/3) and S 3 ∼ (3, 3, 1/3), which are accompanied by a scalar diquark ω ∼ (6, 1, 2/3). The muon receives a large contribution towards its anomalous magnetic moment due to chiralityenhanced effects from leptoquark S 1 that explains a µ data. This same leptoquark S 1 also accommodates for the R D ( * ) anomaly, whereas leptoquark S 3 is responsible to account for the tension observed in the R K ( * ) ratio. Furthermore, both S 1 and S 3 leptoquarks in association with the diquark ω participate in generating masses for the neutrinos at the two-loop order. A detailed analysis is carried out in this work, which shows strong correlations among various flavor violating processes, including neutrino oscillation parameters. In addition to exploring different regions in the parameter space of the theory, we have demonstrated the feasibility of this framework by providing benchmark points. These benchmark points successfully accommodate all three anomalies and naturally incorporate correct neutrino masses and mixings while evading a number of experimental constraints from lepton flavor violation and flavor changing processes, as well as direct searches for leptoquakrs and diquarks at colliders. The lepton flavor violating rare decays of tau lepton τ → µγ, τ → µφ, and τ → µη are all predicted to be right below the current experimental upper bound. Other lepton flavor violating meson decays B s → τ µ and B → Kτ µ are expected to lie around one order below the present experimental limit as well. Hence, this model is very predictive and has the potential to be tested in near future by LHCb and Belle-II. Besides, the presence of TeV scale leptoquarks can lead the way to probe this model at the LHC in near future.