A flavor-inspired radiative neutrino mass model

One of the most important discoveries in particle physics is the observation of nonzero neutrino masses, which dictates that the Standard Model (SM) is incomplete. Moreover, several pieces of evidence of lepton flavor universality violation (LFUV), gathered in the last few years, hint toward physics beyond the SM. TeV-scale scalar leptoquarks are the leading candidates for explaining these flavor anomalies in semileptonic charged and neutral current B-decays, the muon, and the electron magnetic dipole moments that can also participate in neutrino mass generation. In this work, we hypothesize that neutrino masses and LFUV have a common new physics origin and propose a new two--loop neutrino mass model that has the potential to resolve some of these flavor anomalies via leptoquarks and offers rich phenomenology. After deriving the neutrino mass formula for this newly-proposed model, we perform a detailed numerical analysis focusing on neutrino and charged lepton flavor violation phenomenology, where the latter provides stringent constraints on the Yukawa couplings and leptoquark masses. Finally, present and future bounds on the model's parameter space are scrutinized with exemplified benchmark scenarios.


Introduction
The observation of neutrino oscillations was the first direct hint that the Standard Model of particle physics is imperfect and must be extended. Lepton flavor universality (LFU), a solid prediction of the SM, can be easily violated in the beyond SM (BSM) models, where the particles preferentially couple to certain generations of leptons. In the last several years, indications of LFU violation (LFUV) have been observed in both b → s and b → c ν processes. Observables associated with these transitions are the well-known R K ( * ) and R D ( * ) ratios, respectively. LFU in the SM predicts the former ratio to be unity with uncertainties less than 1%. A deficit in this neutral-current transition has been observed consistently over the years in several experiments, and LHCb recently updated their measurements [1] that increased the significance of the deviation. Moreover, with yet another observed deviation in Br B 0 s → µ + µ − [2][3][4][5], the combined significance of the deviation is uplifted to 4.7σ. On the other hand, the R D ( * ) ratio differs from unity due to the substantial mass difference between tauon and muon. An enhancement of this charged-current transition is reported by several experimental measurements, which combinedly leads to approximately 3σ deviation from the SM value [6,7].
Besides, there has been a longstanding tension between the theoretical prediction of the anomalous magnetic dipole moment (AMDM) of the muon (g − 2) µ and the value measured at the BNL E821 experiment [8]. The FNAL E989 experiment [9] has recently announced its result, which has a smaller uncertainty and is fully compatible with the previous best measurement. Together, these two experiments show a remarkably large deviation with a significance of 4.2σ with respect to the theory prediction [10]. Various new physics models are proposed to explain the observed significant departure. For a most recent review see Ref. [11]. The SM prediction given in Ref. [10] is based on the estimate of the leading-order hadronic vacuum polarization contribution, evaluated from a data-driven approach. On the other hand, if recent lattice computations [12][13][14] are considered, then the tension reduces to 1.5σ from 4.2σ. However, if these new lattice results hold, they point towards a large ∼ 4.2σ discrepancy with the low-energy e + e − → hadrons cross-section data with respect to SM predictions [15][16][17][18].
On top of that, the electron AMDM (g − 2) e is also measured in the experiments with an unprecedented level of accuracy. Recently, improved measurement [19] of the finestructure constant utilizing Caesium atom shows a −2.4σ deviation in comparison with the direct experimental measurement [20]. Lately, these anomalies in the lepton AMDMs have gained a lot of attention in the theory community; for simultaneous explains of the muon and the electron AMDMs in various BSM frameworks, see e.g. Refs. . It is noteworthy to point out that a more recent measurement of fine-structure constant utilizing Rubidium atom [73] shows somewhat consistent with the direct measurement of a e [20]. This new result [73] finds ∆a e = +1.6σ, indicating a ∼ 5σ disagreement between these two experiments ( [19] and [73]). Therefore, the electron g − 2 situation requires clarification from future experiments.
All these flavor anomalies mentioned above are strongly pointing toward physics beyond the SM. Interestingly, the prime candidates to solve these flavor anomalies are leptoquarks (LQs), i.e., hypothetical particles that combine the properties of leptons and quarks (for a recent review on LQs, see Ref. [74]). The existence of LQs are highly motivated since particles of this type are naturally predicted by Grand Unified Theories. Explanation of flavor anomalies  requires these particles to have masses of order TeV. Remarkably, these TeV scale scalar leptoquarks (SLQs) can also participate in neutrino mass generation.
In this work, we hypothesize that neutrino masses and LFUV have a common new physics origin. Motivated by this unified framework, we propose a new radiative neutrino mass generation model where scalar leptoquarks, at the leading order, induce tiny neutrino masses as two-loop quantum corrections. If these LQs reside close to the TeV scale, in addition to incorporating neutrino oscillation data, the proposed model has the potential to address the flavor anomalies mentioned above. However, addressing flavor anomalies demands some of the Yukawa couplings to be of order unity. Therefore, with the TeV-scale LQs, charged lepton flavor violating (cLFV) processes are inevitable, which are also clear signals of new physics. In what follows, we first provide the details of our model and derive the neutrino mass matrix in a general gauge. From the derived neutrino mass formula, we carry out a comprehensive phenomenological study of the neutrino sector as well as cLFV, which provides the most stringent constraints on the model parameters. Specifically, we investigate a few minimal benchmark scenarios with a limited number of Yukawa parameters without assuming any strong hierarchy among them. We scrutinize these textures for their ability to satisfy neutrino observables and assess cLFV processes with a detailed numerical study using Markov chain Monte Carlo analysis. Finally, we illustrate how (g − 2) µ , the most prominent flavor anomalies, can be addressed while satisfying all LFV constraints and neutrino oscillation data. This paper is organized in this way: our newly-proposed model of neutrino mass is introduced in Section 2, and the detailed derivation of the neutrino mass formula is given in Section 3. In Sec. 4, we work out in details all charged lepton flavor violating processes that occur in this model and present the results in Sec. 5. Finally, we give our conclusion in Section 6.

Proposed model
The new neutrino mass model proposed in this work consists of the following three BSM scalar multiplets: Numbers in parentheses stand for quantum numbers of each field under SU The SM Higgs is denoted as H(1, 2, 1/2) = H + , H 0 T . This scalar field will get a nonzero vacuum expectation value (vev) v ≡ H = 174 GeV during the spontaneous breaking of the electroweak (EW) symmetry. The fermion sector does not change. It contains the same particle content as in the SM Figure 1: New physics operator leading to non-zero neutrino masses in our proposed model. Here f i represents SM fermions. From the left-most (right-most) vertex, it is clear that R 2 (S 1 ) carries F = 0 (F = −2), where F = 3B + L is known as the fermion number.
where i indicates generation index and ψ c ≡ Cψ R T denotes the charge conjugate of the right-handed field. Among three new scalars introduced, only R and S can be considered as LQs. They interact with the SM fermions through the following Yukawa interactions In order to avoid unnecessary cluttered notation, we have used "·" to denote an SU (2) contraction, e.g., L · Q ≡ L a Q b ab with being the antisymmetric tensor ( 12 = − 12 = 1) and a, b = 1, 2 being SU(2) indices.
All terms in Eq. (2.5) conserve both baryon (B) and lepton (L) numbers. This can be seen, for instance, by setting (B, L) to (1/3, −1) and (−1/3, −1) for R and S fields, respectively. Such assignments forbid diquark terms, QQS † and u c d c S, despite being allowed by the gauge symmetry. Note that it is important to have a globally conserved B, or else a rapid proton decay will take place in our theory. On the contrary, lepton number L is broken, as required to generate non-zero neutrino mass, by the following non-trivial terms in the scalar potential: (2.6) One should have noticed that terms in Eq. (2.6) still conserve the baryon number with ξ field carrying opposite (same) baryon number as that of S (R) field. In the above equation, the two parameters, i.e., λ and µ, can be made real by absorbing their respected phases into scalar fields. The ∆L = 2 effective operators, depicted in Fig. 1, must contain the product of ∆F = 0 and ∆F = 2 couplings, i.e., any combination of y L,R f L,R , so there are four kinds of ∆L = 2 operators that can be generated within this model after integrating out heavy scalar states. It is required that such operators involve gauge bosons, or else they will vanish by SU (2) symmetry. The four operators will contain (HD µ H) multiplied by the following combinations: Note that the SU (2) contraction occurs on fields inside parentheses. In addition, the covariant derivative can also act on other fields. All but the operator (iv) will lead to neutrino masses at two-loop level.
Scalar terms in Eq. (2.6) will cause mixing among ξ 1/3 −S 1/3 , ξ 2/3 −R 2/3 , and ξ 5/3 −R 5/3 components. The mass matrices of leptoquarks relevant for neutrino mass generation are where m R and m ξ are the bare masses of R and ξ, respectively. These two mass matrices can be diagonalized by performing the following rotations where Here c x , s x stand for cos x, sin x. In terms of scalar mass parameters, the two mixing angles are given by Furthermore, the mass eigenvalues of χ 2/3 1,2 and χ 1/3 1,2 are found to be Note that M 2 1,2 and M 2 3,4 can be the larger or the smaller of the two mass eigenvalues. They are defined such that (2.14) In terms of mass eigenvalues, we come out with an alternative way of writing Eq. (2.11), namely from which both λ and µ can be written in terms of mass eigenvalues Having rotated the scalars into their mass eigenstates, we now do the same for Yukawa interactions of Eq. (2.5). Without loss of generality, we can define all couplings in Eq. (2.5) in the charged lepton mass diagonal basis. If we assume further that the up-type quarks be diagonal as well, Eq. (2.5) becomes where V is the Cabibbo-Kobayashi-Maskawa mixing matrix. Note that the results are not affected by changing the up-type to the down-type diagonal basis. Since two choices of the Yukawa coupling textures, namely, "up-type" and "down-type" mass-diagonal basis are widely used in the literature, in the following text, we provide the neutrino mass formula in both these scenarios. As one can see from Fig. 1, the neutrino mass generation requires interaction between LQs and W boson, originating from the SU (2) covariant derivatives where g, g are the SU (2) and U (1) Y gauge couplings and σ a are Pauli matrices. Based from Eq. (2.19), the ξ 2/3 -ξ 1/3 -W vertex can be derived from the triplet kinetic term, that is, After rotating the corresponding fields to their mass eigenstates, we obtain In this model, we work in the general R ξ gauge, so we need to know the LQ interactions with the Goldstone boson. By using Eqs. (2.6) and (2.15), the LQs-Goldstone interactions are found to be will be three subgroups contributing to neutrino masses. Each contribution is presented in Fig. 2. Since neutrino masses are Majorana in nature, in addition to the diagrams shown, there is another set of diagrams with internal particles replaced by their charge conjugates. The sum of the two sets of diagrams will result in the neutrino mass matrix being symmetric. As mentioned before, in evaluating the neutrino mass diagrams, we work in the R ξ gauge. Therefore, the neutrino mass diagrams will contain gauge parameter ξ dependent terms (not to be confused with ξ multiplet). Such terms will later disappear after we sum over all diagrams. The resulted neutrino mass matrix in the up-type quark diagonal basis can be written as Here the factor of 3 accounts for the exchange of color states inside the loops, whereas D u and D are the normalized mass matrices of up-type quarks and charged leptons, re- Similarly, in the basis where down-type quark mass matrix is diagonal, we have One should note that Eqs.
ij denoting the dimensionless loop function for the n-th diagram, that is, The cancellation of terms containing the gauge parameter ξ can be inferred directly from Eqs. (3.5)- (3.12). To see how this cancellation takes place, it is desirable to express k · (2q + k), which appears in all W -mediated diagrams, as Terms inside parentheses will cancel the LQ propagators, particularly those appearing in diagrams 1, 3, and 5. Thus, they will vanish by the orthogonality of LQ mixing matrices. The remaining terms, which are proportional to M 2 b −M 2 a+2 , will make such loop terms have the same coefficients but opposite signs with the corresponding Goldstone loop integrals, allowing the cancellation of ξ-dependent terms. For diagram 7, due to momentum switch between χ 1/3 and χ 2/3 (see Fig. 2 The cancellation of ξ-dependent terms in this case too can be foreseen right away. The gauge parameter cancellation indicates further that all two-loop diagrams presented in Fig. 2 are the complete set of diagrams generating neutrino masses at the lowest order.
We are, then, left with gauge-independent terms. It is straightforward to evaluate the integrals, from which we get (3.14) where only terms relevant to neutrino masses are kept. In those loop integral expressions, we have introduced a parameter All loop integrals are finite, and thus can be calculated numerically. In addition, the contributions of light fermion masses are negligible. Therefore, they can be simply omitted from the integrals, which is demonstrated in Fig. 3.

Lepton flavor violation
All couplings presented in Eq. (2.5) naturally generate lepton-flavor violating (LFV) processes, which are strongly constrained. In this section, we will use those constraints to In this notation, we define as the charged leptons with Q = −1 and subscripts L, R on λ's indicate the chirality of such fields. It is then straightforward to see that there are three types of λ's within this model. The mapping of λ L,R into couplings given in Eq. (2.5) is shown in Table I.
Due to flavor-violating nature of Eq. (2.5), lepton flavor violating processes in general are expected to occur within this model. The first process we consider is i → k +γ * transition, whose effective Lagrangian where F αβ = ∂ α A β −∂ β A α is the electromagnetic field tensor and q = p i −p k is the photon momentum transfer. At the lowest order, this kind of processes arises through penguin-type diagrams exchanging LQs. It is worth noting that, owing to the Ward-Takahashi identity, only dipole terms survive in a lepton decay into an on-shell photon. Its decay width is given by [180] Γ where α em = e 2 /4π is the electromagnetic fine-structure constant, m i is the decaying lepton mass, and In the above equation, 3 is the color factor and x j = m 2 q j /M 2 φ . The summation is performed over all possible couplings and leptoquark fields, as given in Table I. Quantities Q q and Q φ are the corresponding quark and LQ electric charges, which obey Q + Q φ = Q q with Q = −1. Functions F 1 (x) and F 3 (x) are evaluated from diagrams emitting a photon from the quark line, whereas F 2 (x) and F 4 (x) from diagrams emitting a photon from the LQ line. They all are given by Due to the possibility of simultaneous existence of both λ L and λ R in this model, see Table I, we can have a chirality-enhanced process, especially when the top quark is inside the loop. This will lead to severe constraints on the Yukawa couplings. Current and future rates of this kind of processes are presented in Table II.

Lepton 3-body decay
This kind of processes also occurs at loop level, consisting of photon-and Z-penguin diagrams as well as the box diagrams. For the photon-mediated processes, one just needs to attach the photon leg in Eq. (4.2) with a + l − l pair. Now, the photon is off shell, both A 2L,R and A 1L,R contribute to the photon-induced effective Lagrangian, written as [185] L γ−penguin It is also straightforward to evaluate A 1L , A 1R , which are given by with In addition to the aforementioned photonic diagrams, one can also have Z-penguin interactions. They are given by In deriving Z L,R , we have neglected terms proportional to m i . Here g (f ) L,R = T 3f L,R − Q f sin 2 θ W , with θ W being the weak mixing angle, and Similarly, for the box diagrams, we have Note that we do not list box operators in the form of (S ∓P )×(S ±P ), with S, P indicating scalar and pseudoscalar bilinears. This is because they can always be Fierz reordered into the form of (V ∓ A) × (V ± A) operator, which is already included. The corresponding Wilson's coefficients are found to be One can see that B 2L has two terms, but they do not come from the same set of couplings. The first term, coming from momenta of internal quarks, similar to B 1L , has (V − A) × (V + A) type. The second one comes through the internal quark chirality flip, which is then Fierz reordered. That explains why it picks the factor of −1/2. The loop functions b 1 and b 2 are determined to be . (4.14) The factor of −1/2 comes from k µ k ν = 1 2 g µν k 2 , which is later Wick rotated. It is straightforward to evaluate these integrals, yielding Combining all interactions mentioned before, we can write the most general +g 5 (¯ i γ α P R k )(¯ l γ α P L l ) + g 6 (¯ i γ α P L k )(¯ l γ α P R l ) + h.c., ∼ 10 −9 [184] τ + → µ + e − e − 1.5 × 10 −8 [189] ∼ 10 −9 [184] Table III: Current experimental bounds on the BR( i → k m n ). Future sensitivities are presented on the last column.
where we have followed the notation of Ref. [185]. The coefficients g 1 , ...g 6 consist of all contributions from photon, Z, and box diagrams, which are given by From here we can calculate the decay width [185,186] (4.20) In the case of − i → + k − n − n decay (i.e., k = n), only box diagrams contribute. The decay width is found to be [186] where we have definedg i = g i | A 1L,R =Z L,R =0 for i = 3, ..., 6. We present current bounds and projected sensitivities of these processes in Table III.

Lepton anomalous magnetic-dipole moment
This quantity arises from the following interaction with q = p 2 − p 1 and ∆a = F (q 2 = 0), evaluated at one-loop penguin diagrams. This gives (4.23)

µ-e conversion in nuclei
For this process, we are interested in the so-called coherent processes, that is, no change in nucleon state during the transition happens. The relevant interactions, therefore, can be written as [190] L eff = − m µμ σ αβ e LSμ P L e qq + g (q) RVμ γ α P R e + g (q) LVμ γ α P L e qγ α q + h.c. The corresponding Wilson's coefficients are given by µ are given in [190]. The effective couplingsg  Table IV.

Results
In the previous section, we derived all charged lepton flavor violating processes that can be used to constrain the model's parameter space. Since addressing flavor anomalies demands TeV-scale leptoquarks with some of the couplings being of order one, bounds from lepton flavor violating processes provide the most stringent constraints on the Yukawa couplings, which we explore in this section in great details.

Case studies
The neutrino mass formula given in Eq. (3.1) (or Eq. (3.3)) consists of four different Yukawa couplings y L,R and f L,R , which are a priori arbitrary 3×3 matrices. The parameter space is quite broad; therefore, we choose a few specific benchmark scenarios and perform a detailed numerical analysis. Since the terms in the second and the third lines in the neutrino mass formula Eq. (3.1) (or Eq. (3.3)) are proportional to m τ /m t , for Yukawa couplings of a similar order, these terms can be completely neglected. This is why, for our numerical study, we stick to the simplified scenario where f R and y R also provide sub-leading contributions to LFV unless otherwise explicitly mentioned. To further reduce the parameters, we assume vanishing Yukawa couplings with the first generation quarks, i.e., y L 1i , f L 1i = 0. Among the few predictive cases that we consider, in the following, we first discuss the most minimal scenario consisting of six non-zero Yukawa parameters y L 3j , f L 3j = 0 that provides an excellent fit to the neutrino oscillation data. As will be discussed later in the text, further parameter space reduction fails to fit neutrino observables with their respective 2σ values. Considering the loop integral behavior discussed above and the suppression of m τ /m t , in the case of no large hierarchy among Yukawa couplings, it is an excellent approximation to keep only the third generation of quarks. Then the neutrino mass matrix formula, in this case, becomes The above formula applies to both the up-quark and down-quark mass diagonal bases. This is because, in the limit we are working, in the up-quark mass diagonal basis Eq. (3.1), the loop integrals are flavor independent. On the other hand, in the down-quark mass diagonal basis Eq. (3.3), the remaining factor is V tb ≈ 1. One can see from this formula that the neutrino mass matrix is reduced into a rank two matrix, whose determinant vanishes. Thus, this specific texture predicts that one of the neutrinos is massless, although both neutrino mass orderings, i.e., normal hierarchy (NH) and inverted hierarchy (IH), can be admitted. The neutrino mass matrix in this form can nicely fit oscillation data. Before diving into numerics, we first demonstrate that the undetermined Yukawa couplings appearing in the above neutrino mass formula can be fully expressed in terms of neutrino observables and as a function of LQ masses and mixing parameters. By following the parametrization described in [195,196] (for alternative parameterizations, see also, Ref. [197,198]), we determine these Yukawa couplings appearing in Eq. (5.1). To do so, the neutrino mass matrix is diagonalized as follows: where U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix and m i are neutrino mass eigenvalues given by for NH and for IH. Now the neutrino mass matrix given in Eq. (5.1) can be re-written as Here, Y a = y L 3i and Y b = f L 3i are row matrices. Utilizing this form, the two unknown Yukawa coupling matrices can be entirely determined by the known values [199] of neutrino observables, SM fermion masses, and as a function of scalar masses and mixings that run through the loops. For NH, we have Similarly, for IH, the solution for Y a,b takes the following forms In both hierarchies, we define r i = (m i /m) 1/2 . This parametrization is sometimes useful to fix the undetermined Yukawa parameters of the theory. In our detailed numerical analysis, we consider not only the benchmark (BM) scenario as mentioned above but also a few variations of it that include reducing as well as extending the number of parameters. Particularly, all case studies we examine are summarized in the following: • Texture given in Eq. (5.1) with f L 3j , y L 3j = 0. We study both NH and IH, which we label as NH-I and IH-I, respectively. Both cases provide a good fit to neutrino data.
• A more minimal variation of the scenario mentioned above is to choose at least one of f L 3j = 0 or y L 3j = 0, leading to vanishing (M ν ) jj . Coupled with the fact that the lightest neutrino is massless, this restriction clearly does not work for NH. Interestingly, for IH, the case with f L 32 = 0 or y L 32 = 0 can still be fitted within 3σ experimental values of the neutrino observables. We demonstrate this by choosing f L 32 = 0 and label it as IH-II.
• In the cases mentioned above, with all zero entries in the first and the second rows, the (31)-entries of both coupling matrices are required to be comparable with the other entries to provide a good fit, which subsequently leads to large µ → eγ. Consequently, these cases demand large LQ masses to be consistent with the non-observation of LVF.
In search for a minimal texture that is also compatible with ∼ O(1) TeV LQs, we explore a scenario with f L 31 = 0 but introduce nonzero couplings f L 2j , y L 2j for at least one j. Now, although (M ν ) 11 = 0, the determinant of the neutrino mass matrix is no longer zero, so a viable neutrino fit, which is compatible with NH, can be obtained. (A vanishing (11)-element of neutrino mass matrix cannot be realized in IH case.) For demonstration purpose, we choose to have nonzero f L 23 and y L 23 , while other y L 2j , f L 2j are simply set to zero. The two working benchmarks are labeled as NH-II (up-diagonal basis) and NH-III (down-diagonal basis).

Numerical analysis
Our numerical study is based on χ 2 analysis, and the χ 2 -function is defined as where σ i represents experimental 1σ uncertainty; T i and E i represent the theoretical prediction and the experimental central value for the i-th observable, respectively. In the above equation, i is summed over five observables: two neutrino mass squared differences and three mixing angles. For the simplicity of our work, we consider all parameters to be real; hence, we do not attempt to fit the CP-violating Dirac phase in the neutrino sector,  Table V: Neutrino oscillation parameters taken from Ref. [199]. Here, ∆m 2 31 > 0 for NH and ∆m 2 32 < 0 for IH. Here 'bfv' represents best fit values obtained from global fit [199].
which can be trivially done by turning on phases of these couplings. Neutrino oscillation data used in our fit are summarized in Table V. Once a good fit to data is obtained from χ 2 analysis, we perform a Markov chain Monte Carlo (MCMC) analysis to explore the parameter space (consistent with neutrino observables) and inspect lepton flavor violation for which we varied the non-zero Yukawa couplings and LQ masses in the ranges [−1, 1] and [1,100] TeV, respectively. Sample fits obtained from our numerical procedure that is consistent with neutrino observables are presented below for each of the cases listed above (here, we have defined m 0 = a 0m ):    From our detailed numerical scan over the parameter space using MCMC analysis, we obtain interrelationships among various observables: correlations between neutrino mixing parameters, bounds on LQ masses from LFV processes, and correlations among different LVF processes are presented in Figs. 4-9. For both the NH and IH cases, a global fit to neutrino oscillation data has two local minima for the mixing angle θ 23 , the one with θ 23 > 45 • being the lower one [199] (for those without SuperKamiokande data). Whereas for NH, these two minima are almost identical (in the sense of ∆χ 2 measure), however, they significantly differ for IH, and θ 23 > 45 • case is highly preferred to θ 23 < 45 • . This feature is clearly visible in the upper left panel in Fig. 4 for the texture IH-I.
As discussed above, the most minimal Yukawa texture in this theory, which is still consistent with oscillation data, corresponds to IH-II. Due to its minimality, this scenario fails to reproduce neutrino observables within the experimental 2σ range; however, it can be fitted within 3σ values, as can be seen from the third column in Table VI. For this specific texture, a tension exists to simultaneously fit θ 12 and θ 23 close to their central values. This attribute is demonstrated in the upper right panel in Fig. 4.
Moreover, for NH-I and NH-III, both θ 23 > 45 • and θ 23 < 45 • are equally preferred (see middle left and lower panels in Fig. 4, respectively), whereas, for the texture NH-II, MCMC analysis returns solutions only for θ 23 > 45 • as depicted in the middle right panel in Fig. 4.
To obtain a good fit to data, for textures with IH-I, IH-II, and NH-I, the (31)-entries in f L , y L are required to be sizable and are of similar order compared to other non-zero entries, as can be seen from fits Eqs. (5.13)- (5.17). Due to this requirement, the LQ masses  must be much above the TeV scale to satisfy the stringent LFV processes; the most relevant process is the µ → eγ. The plots of this process, as a function of LQ mass, are presented in Fig. 5; they show that M LQ O(10) TeV must be satisfied. On the contrary, for textures NH-II and NH-III, f L 31 is set to zero, and non-zero (23)-entries are introduced. A successful fit to data requires (23)-entries being dominant, whereas y L 31 is somewhat small, as can be seen from fits Eqs. (5.19)-(5.21). Consequently, the branching ratio of µ → eγ is highly suppressed in the latter two scenarios allowing for TeV-scale LQs. However, LQ mass below a TeV is ruled out, and the lower bound on its mass comes from the most dominating LFV processes µ → eee and µ − e conversion, as depicted in Fig. 6.
Further correlations among most prominent LFV processes, namely µ → eγ, µ → eee, and µ − e conversion are depicted in Figs. 7-9. These plots are made by marginalizing over all relevant parameters (LQ mass and Yukawa couplings) in our MCMC likelihood analysis, as discussed above. Remarkably, from these plots, it can be seen that most of the minimal textures we have exploited in this work will be entirely ruled out by upcoming low-energy experiments searching for LFV for M LQ O(100) TeV.
Finally, we demonstrate how to simultaneously satisfy neutrino observables and muon g − 2, where, currently, (g − 2) µ is the most prominent flavor anomaly that shows 4.2σ deviation from the SM prediction. As explained above, textures IH-I, IH-II, NH-I do not allow TeV scale LQs, therefore, to obtain a viable scenario, we consider an example with NH-II texture.
NH-II consists of zero f L 31 entry, so adding a new coupling f R 32 , for instance, will not induce a new contribution to µ → eγ arising from chirality-enhanced term; it only induces a weaker τ → µγ process. However, the benchmark provided in Eq. (5.19) is still unsuitable for incorporating (g − 2) µ with a TeV-scale LQ. This particular fit has a somewhat small f L 32 element; thus, to reproduce (g − 2) µ , an order unity f R 32 coupling is needed. Once this required size of f R 32 is included, along with the f L 33 coupling present in Eq. (5.19), topquark chirality-enhanced contribution to τ → µγ rate becomes too large and rules out this  We verified that introduction of non-zero f R 32 still provides sub-leading contribution to the neutrino mass, and its effect can be safely neglected. However, it has significant effect on cLFV, which we have also incorporated. The above parameters provide a good fit to the neutrino observables Since top-quark chirality enhancement is required to fit (g−2) µ consistently with sizable f L,R 32 entries, as can be seen from Eq. (5.29) that (33)-entry must be pretty small compared to the rest of the elements to keep τ decays under control. In addition, µ − e conversion in the gold nucleus also lies just below the current bound. Branching ratios of these two leading processes for this fit are found to be BR(τ → µγ) = 1.1 × 10 −8 , BR(µ − e) conv. = 6.1 × 10 −13 . (5.33) On the other hand, the consistency with the recent lattice results that weakens the long-standing discrepancy in (g − 2) µ between experiment and theory can be obtained by reducing the value of f R 32 without affecting the neutrino observables. For example, setting f R 32 = 0.134 instead of 0.18 leads to ∆a µ = 1.95 × 10 −9 in agreement with lattice result [12]. Such a reduced value of this coupling subsequently decreases τ → µγ rate; this new value of f R 32 corresponds to BR(τ → µγ) = 5.8 × 10 −9 , whereas BR(µ − e) conversion as quoted in Eq. (5.33) remains unaltered since this coupling plays no role for this observable.

Non-standard Neutrino Interactions
The LQs R 2 and S 1 couple to neutrinos and quarks (cf. Eq. (2.5)), consequently, chargedcurrent non-standard interactions (NSI) at tree-level can be induced [200][201][202]. Using the effective dimension-6 operators for NSI introduced in Ref. [200], the effective NSI parameters in our model can be written (in the up-quark mass diagonal basis) as, whereŷ L ≡ −V T y L . Note that any nonzero y L dα is in conjunction to Cabibbo rotation and inducesŷ L se leading to strong constraints, for instance, K + → π + νν with Re[ŷ L deŷ L se ] = [−3.7, 8.3] × 10 −4 (M S 1 /TeV) 2 . Thus, NSI induced from S 1 LQ via y L Yukawa coupling is subdominant. Moreover, any Yukawa couplings to electron and muon sector f L uα and y L dα (α = e, µ) are subjected to stringent constraints from the non-resonant dilepton searches [149,155] at the LHC. However, the LHC limits on the LQ Yukawa coupling in the tau sector are weaker and in principle be O(1) leading to τ τ as large as 34.4 % [202], which is within reach of long-baseline neutrino experiments, such as DUNE [203].

Conclusions
Neutrino oscillations were discovered almost 25 years ago, showing that neutrinos have a mass; however, its origin remains unknown. Recently, several pieces of evidence of lepton flavor universality violation strongly indicate physics beyond the SM. Scalar leptoquarks are the prime candidates for resolving all these flavor anomalies. Motivated by this, in this work, we hypothesized that neutrino masses and flavor anomalies have a common new physics origin and proposed a new two-loop neutrino mass model consisting of scalar leptoquarks (3, 1, 1/3) and (3, 2, 7/6) along with a third scalar (3, 3, 2/3). Each of these scalar leptoquarks has the potential to incorporate R D ( * ) , (g − 2) e , and (g − 2) µ anomalies. The scalar leptoquark (3, 2, 7/6) may also address anomalies in the R K ( * ) ratios via new physics interactions with the electron. Since resolution to flavor anomalies requires TeV scale scalar leptoquarks with some of the Yukawa couplings of order unity, the proposed model can be tested in ongoing and future colliders. However, probes of lepton flavor violation in neutrino mass models provide the most efficient way of searching for physics beyond the Standard Model that expands far beyond the reach of colliders such as the LHC. In this work, we have primarily focused on the neutrino phenomenology and examined various minimal textures of the Yukawa coupling matrices that can satisfy the neutrino oscillation data. In particular, we have exploited five benchmark scenarios with a limited number of Yukawa parameters of a similar order, two of which provide an inverted hierarchy for the neutrino masses, and the rest provide a normal hierarchy. Moreover, by performing a detailed numerical procedure, namely, the Markov chain Monte Carlo analysis, we have studied in depth various lepton flavor violating processes and constrained the parameter space of this theory. Our analysis shows that for the minimal Yukawa textures considered in this work, the current low-energy experiments provide stringent constraints on model parameters, and near-future experiments hunting for lepton flavor violating rare processes will rule out these scenarios for leptoquark masses below 100 TeV. Finally, we have presented a case study where neutrino observables and the tension in the muon anomalous magnetic moment, the most prominent flavor anomaly, are incorporated simultaneously for TeV scale leptoquark masses by keeping lepton flavor violations under control, which requires a bit of tuning of the Yukawa parameters.