Marriage between neutrino mass and flavor anomalies

Experimental hints for lepton flavor universality violation in beauty-quark decay both in neutral- and charged-current transitions require an extension of the Standard Model for which scalar leptoquarks (LQs) are the prime candidates. Besides, these same LQs can resolve the long-standing tension in the muon and the recently reported deviation in the electron $g-2$ anomalies. These tantalizing flavor anomalies have discrepancies in the range of $2.5\sigma-4.2\sigma$, indicating that the Standard Model of particle physics may finally be cracking. In this Letter, we propose a resolution to all these anomalies within a unified framework that sheds light on the origin of neutrino mass. In this model, the LQs that address flavor anomalies run through the loops and generate neutrino mass at the two-loop order while satisfying all constraints from collider searches, including those from flavor physics.


I. INTRODUCTION
In the Standard Model (SM), the coupling of the electroweak (EW) gauge bosons to leptons is flavor universal, a property known as lepton flavor universality (LFU). The SM's extension is expected to violate this property, and experimental observations of such processes will be direct evidence of physics beyond the SM (BSM). Even though LHC searches have not found any direct hints for new particles, indirect searches for new physics have become an increasingly important pathway. Intriguing hints for BSM physics have accumulated over the years from the measurements of precision observables in several experiments that strongly suggest LFU violation (LFUV).
Lepton anomalous magnetic moments (AMMs) for the first two generations that are very sensitive to new physics are measured in the experiments with unprecedented accuracy. There is a long-standing tension in the muon AMM (g−2) µ measured in 2006 at Brookhaven [1]. Fermilab's new measurement [2] is in complete agreement with the previous result, and the combined result shows a large +4.2σ discrepancy with the SM prediction [3]. For a recent review of BSM models addressing this anomaly, see Ref. [4]. As for the electron, a recent experiment at Berkeley [5] has measured the fine-structure constant using Caesium atom with extreme precision, which indicates a −2.4σ disagreement with the direct experimental measurement [6]. A large positive deviation of (g − 2) µ but a negative deviation of (g−2) e has caused excitement in the high-energy physics community, and various new physics proposals have been made to solve these discrepancies simultaneously .
Here we point out that a more recent measurement utilizing Rubidium atom at the Kastler Brossel Laboratory in Paris [59] shows somewhat consistent with the * E-mail: julio@brin.go.id † E-mail: shaikh.saad@unibas.ch ‡ E-mail: wtd8kz@virginia.edu direct experimental measurement [6]. Contrary to the 2018 result by Berkeley National Laboratory [5], this new 2020 result [59] finds ∆a e to be positive (+1.6σ), indicating a ∼ 5σ disagreement between these two experiments. Since these determinations by the two groups cannot be reconciled within the stated uncertainties, we, at first, will treat these two results as independent cases, presenting separate fits to demonstrate consistency of our model with each result. Later, however, we also perform a fit to illustrate how our model can be simultaneously consistent with the combined results of [5] and [59], favoring a Standard Model-like value. Furthermore, LFU violating B-meson decays have been persistently observed in a series of experiments.
Among these processes, the most prominent deviation has been observed in neutral-current transitions associated with the ratio: R K ( * ) = Br B → K ( * ) µ + µ − /Br B → K ( * ) e + e − , which has very small theory uncertainties since hadronic uncertainties cancel out in the above ratios, making it extremely sensitive to new physics probes. In the SM, this ratio is predicted to be unity [60][61][62][63][64] due to the LFU property. However, experiments have consistently found it to be somewhat smaller than one [65][66][67][68]; among them, the most precise measurement by LHCb [69] alone quotes a significance of 3.1σ for R K -ratio. 1 On the other hand, charged-current transitions in the B-meson decays indicate an enhancement of the ratio: [82][83][84][85][86] from the SM predictions [87,88] with a combined significance of about 3σ [89]. 1 In addition to the clean observable R K ( * ) , several other observables are in tension with the SM predictions for example angular distributions and branching ratios of several b → sµ + µ − modes [70][71][72][73][74][75][76][77][78][79]. When all these anomalies are combined, a global fit [80,81] to data suggests > 5σ discrepancy. In this Letter, we only focus on the R K ( * ) ratio; the rest of the observables suffer from large hadronic uncertainties, and we do not attempt to alleviate these discrepancies.
Besides all this, first and foremost, it is widely accepted that the SM cannot be the fundamental theory at all energies since it fails to incorporate neutrino oscillations firmly confirmed by various experiments; for reviews see Refs. [90,91]. Motivated by this crucial drawback of the SM, in this Letter, we propose a simultaneous solution to all four of the anomalies mentioned above that are directly intertwined with the neutrino mass generation mechanism. In this unified framework [92], we employ two scalar leptoquarks 2 (LQs), {S 1 (3, 1, 1/3), R 2 (3, 2, 7/6)}, which address the flavor anomalies and neutrino masses simultaneously. In this model, Majorana neutrino mass appears at two-loop order. Our novel proposal is minimal, and the only neutrino mass model in the literature that resolves all of these aforementioned notable anomalies while satisfying low-energy flavor constraints and the LHC limits. 3 Furthermore, the proposed model is fully testable since a resolution of the charged-current B-anomaly requires the relevant LQ to have a mass around 1 TeV (the neutralcurrent and g − 2 anomalies can be explained with relatively heavier LQ).
This paper is organized in the following way: we introduce the model in Sec. II and discuss how anomalies are incorporated within this setup in Secs. III and IV. After presenting our numerical results in Sec. V, finally we conclude in Sec. VI.

II. MODEL
In addition to the SM particle content, the proposed model consists of two scalar LQs (SLQs), S 1 (3, 1, +1/3), R 2 (3, 2, 7/6), and another BSM multiplet ξ 3 (3, 3, 2/3). We denote their component fields by The couplings of LQs with the SM fermions take the following form: where Q and L stand for, respectively, left-handed quark and lepton doublets, u R and R are right-handed up-type quark and lepton, and i, j = 1 − 3 are family indices. We also use "·" to denote SU (2) contraction so that L · Q ≡ L ρ Q σ ρσ with being the Levi-Civita tensor and ρ, σ = being SU (2) indices. Since explanations of flavor anomalies require that these LQs couple to quark-lepton bilinears, as usual, we turn off the diquark coupling of S 1 LQ, which guarantees baryon number conservation. Now, expanding these Yukawa terms, we obtain where, V represents the CKM mixing matrix, and we have chosen to work in the so-called "up-quark mass diagonal basis". The two LQs and ξ can also interact with Higgs doublet H in the scalar potential V sc ⊃ λS † 1 H T ξ † H + µR † 2 ξH. The simultaneous presence of the Yukawa and these scalar interactions breaks the lepton number by two units, as required for generating non-zero neutrino mass; one of the diagrams is shown in Fig. 1. If the Yukawa interactions are turned off, the leptoquarks are not required to carry any lepton number; hence the scalar potential will remain lepton-number conserving. After the EW symmetry breaking, these terms induce mixing in the LQ-ξ sector, yielding six physical exotic scalars, i.e., χ 1/3 a , χ 2/3 a , and χ 5/3 a with a = 1, 2. Since explaining flavor anomalies requires some O(1) Yukawa couplings, we are compelled, by tiny neutrino masses, to have small mixing angles. Therefore, it is a good approximation to identify gauge eigenstates R 2/3,5/3 and S 1/3 as their physical eigenstates.
In the chosen basis, the neutrino mass formula takes the following form: where m 0 = 3g 2 m t / √ 2(16π 2 ) 2 with g being the SU (2) gauge coupling constant and m t being the top-quark mass. The normalized mass matrices of up-type quarks and charged leptons are D u = diag. mu mt , mc mt , 1 and D = diag. me mτ , mµ mτ , 1 . The loop integral I 0 in the asymptotic limit (i.e., when all fermion masses are zero) is given by Here θ and φ denote the mixing angles of LQ with charges Process Constraints µ → eγ [198] |f . The second (first) row in the "Constraints" column shows the bound with (without) chiral enhancement. Constraints on the S1 LQ couplings y L,R with no chiral enhancement is weaker by a factor of 3, whereas chirally enhanced contribution has the factor with C replaced by C = 7 + 4 log is the loop function, derived by assuming M b /M a+2 = 1 + r ba with r ba 1. Note that M 1,2 are masses for 2/3 charged LQs and M 3,4 are for 1/3 charged LQ. For more detailed discussion, we refer the reader to Ref. [92]. Figure 1. A typical two-loop diagram leading to non-zero neutrino mass. For full set of diagrams, see Ref. [92].

III. CORRELATED OBSERVABLES
This section discusses the relevant correlated observables/constraints in the proposed two-loop neutrino mass model that resolve the four anomalies.
At tree-level both S 1 and R 2 LQs contribute to the b → cτ ν transition which are favorable, whereas they do not lead to a desired contribution to the b → sµµ. Loop-induced contributions to b → sµµ transition face several difficulties and unable to consistently solve R K ( * ) anomalies (see e.g. Ref. [174]). On the contrary, if R 2 LQ couples to the electrons, relatively small couplings are feasible to address [154,201] R K ( * ) anomalies since the relevant operator is generated at the tree-level. Moreover, each of these LQs can also simultaneously explain both the (g − 2) µ and (g − 2) e anomalies via the chirallyenhanced top-quark and charm-quark contributions, respectively [33].
To be specific,f R se,be (for brevity, anomaly requires either f R,L tµ or y R,L tµ (f R,L t/ce or y R,L t/ce ) couplings to be sizable. Furthermore, with TeV scale LQ (LQ mass below TeV is highly constrained by LHC data), Consequently, the NP contribution to resolving these four anomalies while incorporating neutrino oscillation data is highly non-trivial and leads the way to various flavor violating processes that are severely constrained by the experimental data. One of the most important constraints is the lepton flavor violating (LVF) radiative decay of charged leptons i → j γ [202]. For instance, the µ → eγ process precludes R 2 LQ to accommodate ∆a µ while resolving R K ( * ) . Furthermore, to avoid strong constraints from chirally enhanced i → j γ via the CKM rotation in the up-sector, we have chosen to work in the up-quark mass diagonal basis of Eq. (2). Constraints on the Yukawa couplings owing to dominant LFV processes as a function of LQ mass for both non-chiral and chirally enhanced diagrams are shown in Table I together with constraint from coherent µ − e conversion in nuclei.
Similarly, R 2 LQ in the up-quark diagonal basis leads to various kaon decay processes [203][204][205] for which numerous stringent bounds are listed in Table II. Moreover, the model requires f R bτ ∼ O(1), which modifies Z-boson decay to fermion pairs through one-loop radiative corrections. Then the LEP data [206] on the effective coupling leads to |f R bτ | ≤ 1.21 [207] for LQ mass of 1.0 TeV.
Non-standard neutrino interactions: Both R 2 and S 1 LQs have couplings to neutrinos and quarks that can induce charged-current NSI at tree-level [208][209][210]. Using the effective dimension-6 operators for NSI introduced in Ref. [208], the effective NSI parameters in our

Process
Constraints [203][204][205]   model are given by The Yukawa couplingŷ L dα is subject to strong constraint from K + → π + νν as shown in Table II. The couplings f L uα (α = e, µ) are constrained by the non-resonant dilepton searches at LHC and the limit on f L ue and f L uµ are 0.19 and 0.16 for 1 TeV LQ mass. Thus, the contribution of y L dα and f L uα which leads to ε αβ is at sub-percent level. 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 % [210], which is within reach of long-baseline neutrino experiments, such as DUNE [211].

IV. CONSTRAINTS FROM DIRECT SEARCHES
At the LHC, S 1 and R 2 LQs can be pair produced [212][213][214] through gg and qq fusion processes or can be singly produced in association with charged leptons via s-and t-channel quark-gluon fusion processes. The single production of LQ becomes only relevant for larger Yukawa couplings to the first and second-generation quarks [215]. Thus, collider bounds from single-production are not so significant compared to QCD-driven LQ pair production for our analysis. In both ATLAS and CMS, there are dedicated searches for the LQ pair production in different modes, LQ † LQ → qq ¯ , qqνν. The collider limits on LQ mass depend on the branching ratios to different modes, as depicted in Fig. 2 (upper-plot). Apart from the LQ-pair production bound, there are bounds on the couplings and mass on the LQ from the high-p T tails of pp → distributions [133,140,169,174,[216][217][218][219][220]. The corresponding bounds are also shown in Fig. 2 (lowerplot). The above direct searches assume that the LQs decay promptly, which is no longer applicable for long-lived states. The colored particles arising from the ξ field can generally be long-lived if their mixing with the LQs is sufficiently small, typically required to correctly incorporate the neutrino mass scale. Colored long-lived particles (CLLPs) hadronize prior to decay, forming so-called R-hadrons [221] that are bound states composed of the CLLP and light quarks, antiquarks, and gluons. Since these particles carry color charge, for TeV scale masses, they are efficiently pair-produced pp → χ a χ a via pure QCD interactions. Once produced, they decay to SM fermions with decay width Γ ∼ y 2 m χ /(16π) [222] corresponding to a vertex displacement of d ∼ 10 −14 /y 2 mm. Therefore, for Yukawa couplings y 10 −6.5 , decay is prompt since the corresponding vertex displacements are smaller than typical LHC detector resolution [212]. On the other hand, for y 10 −6.5 they become longlived, and for y 10 −9 these states are practically stable and decay outside the detector. Long-lived gluinos and squarks have been extensively searched for at the LHC. However, no signal above the expected background is found, which provides lower bounds on the CLLPs. Our scenario is similar to the long-lived sbottom (Q em = 1/3)  Table III. 3σ allowed ranges of the neutrino oscillation parameters are taken from a recent global-fit [225]. 1σ allowed ranges of lepton g − 2 and C ee 9 = C ee 10 are also summarized. The last two columns contain the fit values of the observables for each texture. For (g − 2)e, the number in the parenthesis in the second column corresponds to the measurement of [59], and the corresponding fitted values are given in parentheses in the third and fourth columns for TX-I and TX-II, respectively. and stop (Q em = 2/3) for which the current bounds at 95% confidence level are m ≥ 1250 GeV [223] and m ≥ 1800 GeV [224], respectively.

V. RESULTS AND DISCUSSION
In this section, we present our numerical results and study the correlations among R D ( * ) , R K ( * ) , ∆a µ , and ∆a e anomalies within their 1σ measured values while consistently incorporating neutrino oscillations data. First, we demonstrate the viability of our proposed model by providing two benchmark points 4 given in Eqs. (7) and (11)   (bold) from R 2 LQ. Non-zero elements in black are the additional entries required to explain neutrino oscillation data. The only non-zero entry that induces observable NSI is shown in italics.

TX-II:
With m 0 I 0 = 3.33 × 10 −6 GeV, M R2 = 1.0 TeV and M S1 = 1.5 TeV: This flavor structure, on the other hand, explains , and (g − 2) e ( ) by making use of R 2 LQ, and (g − 2) µ ( † ) from S 1 LQ. For each of these textures, fitted neutrino observables are presented in Table III, which is in excellent agreement with experimental data. The consistency of these benchmarks fits with the current LHC bounds are discussed later in this section.
Fits presented in Eqs. (7) and (11) incorporate the (g − 2) e measurement of Ref. [5]. Here we clearly demonstrate how to consistently reproduce the more recent measurement of Ref. [59]. Interestingly, for fit Eq. (7), the coupling y R 21 does not play any role in neutrino mass generation. Hence reducing its size by a factor of ∼ 1/2 and flipping its sign (to be precise, the new value is y R 21 = 0.17) reproduces the result of Ref. [59] without affecting any other observable. On the other hand, for fit Eq. (11) to replicate the measurement of Ref. [59], f R 21 needs to be reduced while keeping its product with f R 31 same. Moreover, the couplings f R 21,31 play an insignificant role in neutrino mass generation as it is always accompanied by electron mass (i.e., suppressed by a factor of m e /m t , see Eq. (3)). For instance, f R 21 = 0.0093 and f R 31 = 0.252 correctly reproduce the fit values illustrated in Table III as TX-II, while reproducing the (g−2) e result of Ref. [59].
It is worth noting that our model can, in fact, be simultaneously consistent with both values of (g − 2) e inferred from Cs [5] and Rb [59] experiments. To illustrate this, we need to combine the two uncorrelated values a e,Cs ± σ Cs and a e,Rb ± σ Rb of the same physical quantity a e to get the best estimate a comb. e ± σ comb. via [228] a comb.
where k = Cs, Rb. For the reader's convenience, here, we collect those values of a e inferred from the measurements of the fine-structure constant: a e,Cs = 0.00115965218161(23) [5] and a e,Rb = 0.001159652180252(95) [59]. Using Eq. (9), we obtain a comb. e = 0.00115965218045 (9). Comparing with the direct measurement, i.e., a exp.  is heavily weighted by a Rb e in Eq. 9). Now, for TX-I with NH, one can turn off Yukawa coupling y 21 R without changing the corresponding fit (see discussion above). This choice gives no new physics contribution to electron (g − 2), i.e., a special case corresponding to a pure SM-like scenario with ∆a e = 0. Furthermore, the combined value given in Eq. (10) can also be trivially reproduced by setting y 21 R = 0.095 without affecting any other observable.
On the other hand, for IH with texture TX-II, setting f L 21 = 0 would lead to the texture with (M ν ) 11 = (M ν ) 22 = 0, forbidden by the oscillation data. One can however, lower the Yukawa coupling f R 21 while keeping its product with f R 31 fixed to obtain C ee 9 = C ee 10 within 1σ. Moreover, f R 21 0.0067 is required to satisfy the bound |f R 31 f R * 33 | < 0.32 (M R2 /TeV) 2 obtained from τ → eγ. This choice of f R 21 leads to ∆a e 4.32×10 −13 , consistent with combined measurement given in Eq. (10).
Next we show that a slight deviation from the TX-II opens up parameter space allowing even a smaller value of ∆a e . We demonstrate the viability by the following modified texture: The above benchmark, which is consistent with all flavor violating constrains, reproduces the central value of the neutrino observables NuFit5.1 [225] and returns same values for the rest of the observables as in TX-II, except for (g − 2) e , which is found to be ∆a e = 7.7 × 10 −14 . This number is also fully consistent with the combined ∆a e given in Eq. (10). Next, we study the correlations among numerous observables associated with the benchmark scenarios, TX-I and TX-II. For TX-I, the allowed parameter space to explain R D ( * ) and R K ( * ) at 1σ (green shaded) and 2σ (yellow shaded) CL in the relevant Yukawa coupling planes are illustrated in Fig. 3 by fixing R 2 (S 1 ) mass at 1.2 (1.5) TeV. In computing R D ( * ) and R K ( * ) observables, we utilized Flavio package [63]. By fixing y R tµ = 0.0025 (cf. Eq. (7)), the allowed range to incorporate ∆a µ at 1σ is shown in orange band that corresponds to y R tµ y L tµ ∼ 10 −3 .  We also show the exclusion regions obtained from resonant (ttµμ) and non-resonant (pp → τ τ, ee) production of LQs, along with flavor constraints such as µ → eγ for a fixed f R bµ , and kaon decays. Here the remaining anomaly ∆a e is determined by the product y R ce y L ce , and the corresponding fitted value is listed in Table III obtained from the fit Eq. (7).
Similarly, Fig. 4 illustrates the TX-II scenario where R 2 LQ addresses R D ( * ) , R K ( * ) , and ∆a e anomalies. Green and yellow bands correspond to 1σ and 2σ allowed ranges for R D ( * ) and R K ( * ) , respectively, whereas orange band shows 1σ preferred value for ∆a e . In this plot, M R2 = 1 TeV and f L ce = 0.3 (cf. Eq. (11)) are kept fixed. The exclusion regions are obtained from resonant (ttττ , cceē) and non-resonant (pp → τ τ, ee) production of LQs, along with flavor constraints such as τ → eγ for a fixed f R bτ and f R be , kaon decays, and Z-boson decay Z → τ τ . ∆a µ , which is not depicted in the figure, is determined by the product y R tµ y L tµ , and the corresponding benchmark value is shown in Table III   As can be seen from these plots, the valid parameter space addressing all the anomalies is very limited due to various collider and flavor violating constraints, and upcoming experiments searching for LFV [198,229] have the potential to rule out this scenario. On the other hand, a resolution to flavor anomalies, particularly the R D ( * ) ratios, requires LQ masses not higher than O(1) TeV, making this model fully testable at future highluminosity LHC.
LHC bounds on fits TX-I and TX-II: For TX-I with f L uτ = 0 (contributing to only NSI), R 2 LQ primarily decays to jjeē and cceē [230] with the same Yukawa coupling f R se leading to branching ratio β = 0.5. Here we fix R 2 mass at 1.5 TeV (corresponding to branching ratio β 0.74 to jjee). Moreover, once f L uτ ∼ O(1) is included, new modes open up, and R 2 primarily decays to jjττ and jjνν [231]. There are no dedicated searches for jjττ and the limits from jjνν is 1.0 TeV for β = 1. On the other hand, for the TX-II, R 2 dominantly decays to ttττ [232], bbττ [233,234], jjττ , and jjνν each with the branching ratio β 0.24. Here we fix R 2 mass at 1 TeV, which allows the corresponding branching ratio to be as large as β 0.27 (from ttττ ).
On the other hand, for TX-I, S 1 primarily decays to jjττ , ttµμ [235], and bbνν [231] with the branching ratios (0.26, 0.37, 0.37), respectively. The latter two provide limits on the scalar LQ to be 1.5 (1.2) TeV and 1.1 (0.8) TeV for branching fraction β = 1 (β = 0.5), which is easily satisfied for our choice of S 1 LQ mass of 1.2 TeV. Similarly, for the TX-II S 1 decays to ttµμ 100% of the time with the bound of 1.5 TeV as previously mentioned. For TX-II, we fix S 1 LQ mass at 1.5 TeV such that current collider constraints along with all the flavor constraints are easily satisfied.
Finally, we discuss the collider bounds on the states mostly originating from ξ. For our fits, using the definition of m 0 (∼ 6 · 10 −3 GeV), we find I 0 ∼ 2 · 10 −4 . Moreover, for all physical scalars with masses close to each other (i.e., r ba 1), utilizing Eq. (4)-(5), we obtain I 0 ≈ 2 sin 2θ sin 2φ. Then, assuming mixings of similar order, we obtain θ ∼ φ ∼ 10 −2 . Therefore, states arising primarily from ξ have Yukawa couplings of order 10 −2 with the third generation fermions (couplings with the first and second generations are more suppressed), hence decay promptly. Consequently, collider bounds on these states are identical to the bounds on LQs.

VI. CONCLUSIONS
In this Letter, we have presented for the first time a two-loop radiative neutrino mass model consisting of two scalars LQs S 1 (3, 1, 1/3) and R 2 (3, 2, 7/6) and another BSM scalar ξ(3, 3, 2/3) with close-knit connections with R K ( * ) , R D ( * ) , and (g − 2) µ,e . We have performed a detailed analysis, including all relevant flavor violating and collider constraints, and presented benchmark fits showing complete consistency with anomalies and neutrino oscillation data. Furthermore, resolving these tantalizing flavor anomalies requires LQs to have masses around the TeV scale, enabling us to test this theory at the ongoing and future colliders. Finally, it is worth mentioning that the model sets several exciting avenues for future work; the implications on electric dipole moments and their correlations with the CP-violating phases in the neutrino sector, detailed analysis on LFV decays, and a dedicated collider study.
Acknowledgments.-The work of J.J. was supported in part by the National Research and Innovation Agency of the Republic of Indonesia via Research Support Facility Program.