Common explanation to the $R_{K^{(*)}}$, $R_{D^{(*)}}$ and $\epsilon^\prime/\epsilon$ anomalies in a 3HDM+$\nu_R$ and connections to neutrino physics

Scalar theories can account for the current $R_{D^{(*)}}$ measurements through a vector operator $\bar{c}_L \gamma_{\mu} b_L\,\bar{\tau}_L \gamma^{\mu}\nu_L$ induced at the loop level. Once the vector contribution is considered on top of a subdominant tree-level scalar component, the predicted value of $R_{D^{(*)}}$ falls within the $1\sigma$ region indicated by the experiments. We explicitly demonstrate this claim in the framework of a three Higgs doublet model extended with GeV scale right-handed neutrinos, by matching the anomalous signal for perturbative values of the involved couplings and respecting the bounds from complementary flavour physics measurements. Remarkably, we furthermore show that the proposed framework can be employed to simultaneously explain also the present $R_{K^{(*)}}$ measurement, as well as the deviation in $\epsilon'/\epsilon$ currently being debated in the literature. These results are obtained by considering the contribution of relatively light right-handed neutrinos which are fundamental in mediating the processes behind the anomalous signals. In this way our findings reveal a new possible connection that links the flavour anomalies to the phenomenology of extended Higgs sector and neutrino physics.

These difficulties are grounded in the different origin of these anomalies, which seemingly require new contributions to appear at several scales and put into discussion properties, such as the universality of gauge couplings, thoroughly verified in collider experiments. For instance, since R D ( * ) quantifies the ratio of processes that receive tree-level contributions within the SM, the measured value of this parameter naturally calls for new degrees of freedom below the electroweak scale. On the contrary, as R K ( * ) is purely determined by the loop structure of the SM, we naively expect that deviations in this observable originate at a much larger scale barring the fine-tuning of the involved couplings. Furthermore, the anomalous signals require new sources of lepton flavour universality violation which, in concrete frameworks, are generally tightly constrained by collider and complementary flavour experiments.
If we analyse the SM case, we see that the only source of flavour non-universality is the Higgs sector, where the Yukawa couplings of fermions necessarily reflect the measured mass hierarchy. Consequently, extending the Higgs sector could be a natural way to implement the additional sources of flavour non-universality at the basis of the measured signals.
In the present paper we pursue this possibility, showing that the R K ( * ) and R D ( * ) measurements can simultaneously be addressed with an extended Higgs sector. The result is based on the preliminary investigation in Ref. [29], where we demonstrated how scalar extensions of the SM can generate a sizeable vector operatorc L γ µ b LτL γ µ ν L at the loop level. Explaining the mentioned flavour anomalies in respect of the complementary experimental bounds requires at least two new scalar particles, which must belong to different gauge multiplets, as well as the inclusion of right-handed neutrinos (RHNs). A natural ultraviolet completion of the proposed framework is then provided by a three Higgs doublet model (3HDM), augmented with RHNs, which we analyse in detail in the following.
As we will show, the inclusion of the tree-level scalar operators contributions allows the model to reproduce the anomalous R D ( * ) signal for values of the involved couplings of order unity. We consider all the emerging effective chiral operators, including the ones mediated by relatively light right-handed neutrinos usually omitted in these studies. We explicitly propose the Yukawa coupling texture of the additional Higgs doublets needed to explanation the anomaly. Remarkably, we find that the same texture allows also for the explanation of the R K ( * ) anomaly owing to the presence of Majorana RHNs, used here exclusively to mediate the processes behind these signals 1 . In this light, flavour anomalies could therefore constitute a first indirect collider signature supporting the existence of RHNs. The proposed framework then connects the issue of flavour measurements to other open questions within neutrino physics and cosmology, such as the neutrino mass generation mechanism and the dynamics of baryogenesis via leptogenesis, extending the phenomenological reach of flavour experiments. We highlight that RHNs in the considered mass range can be discovered by the planned SHIP experiment at CERN [41], as well as in collider experiments through the implied lepton number violating effects [42].
A brief part of the present paper is dedicated to the assessment of a further potential flavour anomaly, connected to the amount of direct CP violation measured in K → ππ decays and quantified in the ratio / . In this case, without advocating for the necessity of a new physics contribution, we find that a minimal modification of the proposed Yukawa texture would allow to explain the signal, should it survive further scrutiny by the dedicated community.
The paper is organized as follows. Section 2 is dedicated to the R D ( * ) anomaly: after a brief review of its experimental status, we show how the considered framework explains the anomalous signal. The R K ( * ) anomaly is studied in a similar fashion in Sec. 3. In Sec. 4 we then present the result of our investigation, analyse the compatibility of the two anomalies and detail the potential impact of new determinations of / . Finally, our conclusions are presented in Sec. 5.

The R D ( * ) anomaly
We start by briefly reviewing the present experimental status of the R D ( * ) anomaly, setting out the formalism at the basis of the analysis.

Experimental status and effective Lagrangian
The results of the LHCb experiment [1,2] have underlined the presence of lepton flavour violating dynamics in the charged decays of the B-meson. In more detail, anomalies have been measured in the ratio of branching fractions for = e, µ. Being defined as ratios, R D and R D * are observables of particular importance that test lepton universality in the transition b → cτν τ net of the involved hadronic uncertainties.
If considered along the initial measurements by the BaBar [3,4] and Belle [5][6][7] collaborations, the LHCb observations set [43] R exp D = 0.407 ± 0.039 ± 0.024 , 1 Alternative explanations of the anomalous signals which rely on the presence of RHNs in the involved final states can be found in Refs. [25,26,[38][39][40] highlighting a strong affinity for third generation leptons. The implied departure from universality falls well beyond the SM predictions with a combined significance of about 4σ [43].
The possible presence of new physics in B anomalies is investigated by detailing the low-energy regime of a potential model, through the construction of an effective theory characterized by dimension six operators that preserve color and electric charge. Then, at the b-quark mass scale, the theory is described by the effective Lagrangian [44][45][46][47] L b→c ν where = e, µ, τ . The four effective operators that appear above are given by In this formalism, the dynamics of the b → c ν process is fully encoded in a corresponding set of Wilson coefficients C = C V L , C AL , C SL , C P L , which parametrize the relevant observables through the involved branching fractions The contributions of heavy states to the effective Lagrangian in eq. (4) is computed by matching the full theory to the effective one at the scale where the heavy degrees of freedom are integrated out. The final expression of the Lagrangian is then obtained upon the computation of the RGE evolution down to the scale at which the process is probed. The matching has been performed diagrammatically and we checked that the QCD and QED corrections induced by the RGE evolution have a negligible impact on the couplings in eq. (13).
A general analysis of Wilson coefficients provides a first portray of the framework we seek. For instance, the simplest scalar extensions of the SM yield tree-level contributions to the scalar and pseudoscalar Wilson coefficients C τ SL and C τ P L which, in principle, allow to explain the anomalous B-physics signal. As shown in the left panel of Fig. 1, new contributions to these quantities let the theoretical predictions enter the 68% confidence interval associated to the signal, denoted by the red areas in the plot. Such a solution, however, is invalidated by measurements of the B c lifetime, which severely constrain the pseudoscalar Wilson coefficient due to the mass hierarchy of the SM: The impact of this constraint is represented in both the panels of Fig. 1 by the areas shaded in light and dark gray, which indicate the values of the Wilson coefficients that result in deviations larger than 10% or 30% from the measured B c lifetime, respectively. As we can see, solutions characterized by large values of C τ P L are obviously disfavored.
As originally proposed in Ref. [29], a possible way for scalar extensions to cope with the B c lifetime constraint is to rely on scalar loop contributions to the vector and pseudovector Wilson coefficients. The case is shown in the right panel of Fig. 1, which makes clear that modest values of C τ V L and C τ AL allow to explain the observed anomaly even when the strongest constraint on the B c lifetime is considered. A potential issue with this solution is that Figure 1: Model-independent fit of the anomalous signal as a function of the indicated Wilson coefficients from eq. (4). The 68%, 95% and 99% confidence intervals selected by the joint fit of R D and R D * are marked by the red, orange and yellow areas, respectively. The shaded light (dark) gray areas show instead the current bound from the measured B c lifetime assuming a 10% (30%) maximal allowed deviation.
contributions to C τ V L and C τ AL are generated in scalar theories only at the loop-level. The required magnitude, of about O(10 −1 ), may then impose couplings on the verge of non-perturbativity [29].
We address here this problem by explaining the observed violation of lepton universality through vector and pseudovector contributions on top of a subdominant pseudoscalar component, used to relax the values of the coupling involved in the loop diagrams that generate C τ V L and C τ AL . The interplay between the involved Wilson coefficients is shown in the right panel of Fig. 1, where crescent values of C τ P L allow an excellent fit of the anomaly for lower values of the vector and pseudovector contributions.

The 3HDM contribution to R D ( * )
In order to respect the tight bounds on B → X s γ, as well as further constraints imposed by penguin diagrams, we follow the setup of Ref. [29] and consider three different scalar SU(2) doublets with H 0 being the SM Higgs doublet. We arrange the scalar potential so that the new doublets do not develop vacuum expectation values, and restrict their masses to the ∼ 300 − 350 GeV range to remain within the reach of current collider searches. As mentioned before, we also consider three RH neutrinos ν R with corresponding Majorana mass terms. In regard of this, we require that m b − m c − m τ m ν τ R to prevent the related b decay channel, maintaining however m ν i R << m t to retain sizeable loop contributions. From the neutrino physics point of view, these additional states are therefore sufficiently heavy to decouple from the low-energy dynamics of active SM neutrinos.
The interactions of the new scalar states are detailed in the following Lagrangian where, for sake of minimality, we set Y d 2 = 0.
The Yukawa texture that we investigate is of the form In the above equations, as well as throughout the rest of the paper, the symbol f denotes the coupling of the H 1 doublet with the fields indicated by the subscript. Analogously, we indicate with g the couplings of H 2 . The elements of the new Yukawa matrices rendered in violet regulate the contribution to the pseudoscalar Wilson coefficient sourced by the first diagram in Fig. 2. The elements in teal enter, instead, the vector and psudovector Wilson coefficient through the remaining loop diagrams. Notice that fτ L τ R is involved in both contributions, depending on the considered component of H 1 . Lastly, terms in orange affect exclusively the computation of the R K ( * ) anomaly presented in Sec. 3. For sake of simplicity we assume real couplings and set to zero the remaining elements of the Yukawa matrices, with the understanding that their values are negligible within our effective description. Figure 2: Additional diagrams for the process b → cτν supported by the considered 3HDM. The first diagram sources the tree-level expression for C τ SL and C τ P L given in eqs. (14) and (15). The last two diagrams, instead, yield the contributions to C τ V L and C τ AL reported in eqs. (16) and (17). We denoted with the symbol ν an active SM neutrino.
By integrating out the degrees of freedom above the b-quark mass scale and matching the Lagrangian in eq. (12) to the effective one in eq. (4), we identify the following tree-level contributions to the scalar and pseudoscalar Wilson coefficient: where g w is the coupling constant of SM weak interactions. Notice that gauge invariance imposes fν L τ R ≡ fτ L τ R . Eqs. (14) and (15) make clear the choice of the consider Yukawa pattern, which allows fb L c R and fc L b R to separately source the tree-level contribution. In this way, two independent degrees of freedom regulate C SL and C P L , making it possible to exploit a CKM enhancement to maintain perturbative values of the involved couplings.
The contributions to the vector and pseudovector Wilson coefficient due to the loop diagram in Fig. 2 where we used gτ L ν R = g ν L ν R and indicated with D dd00 the 4-point loop integral: Clearly, the final expressions for the Wilson coefficients are obtained as C τ V L = C τ (1) AL . Before moving to the discussion of the R K ( * ) anomaly, we address below the main experimental bounds that the proposed solution for R D ( * ) faces.

Main experimental constraints
The main experimental bounds that oppose to the proposed solution for the R D ( * ) anomaly are due to measurements of B → X s γ, B + → K + νν, and B 0 → K * 0 νν.
The introduction of extra SU(2) doublets, through their Yukawa interactions, affects the dimension-5 photon and gluon dipole operators and modifies via the effective Lagrangian [48] the SM prediction for B → X s γ [48][49][50]. It is then clear that measurements of this quantity limit possible new physics contributions resulting from eq. (13). Employing comparable couplings for both the extra doublets to induce an enhancement of the box contribution to R D ( * ) of about a factor of 4, in particular, generates unacceptable large contributions to B → X s γ from the diagrams in Fig. 3. This is manifest even in the limit where one of the two scalar doublets does not participate in the dynamics of the anomaly (2HDM limit 2 ), resulting in the bound on new contribution to C τ V L from B → X s γ plotted in Fig. 4.
The bound therefore provides a clear indication in favour of our mechanism, which relies on separate gauge multiplets with complementary roles that generate the required amount of lepton flavour universality violation without inducing large corrections to the photon and Z vertices. We remark that, within the full model, it is in principle possible to adjust the value of Y d 1 3,3 to cancel the new large corrections to b → sγ without suppressing, at the same time, the box diagrams controlling the anomaly. However, we do not pursue such a possibility because of the amount of fine tuning implied. As we have seen, measurements of B → X s γ yield important constraints on the Yukawa couplings of the new Higgs doublet within the quark sector. In order to investigate similar constraints that potentially target the leptonic sector, we turn now our attention to b → sν LνL . Our 3HDM indeed produces a similar enhancement in the related processes via dimension-6 four-fermion operators encoded in the effective Lagrangian The main contribution is due to the diagram in Fig 5, which leads to a new flavor violating deformation of C τ L given by: Notice that the absence of RHNs in the final states prevents f ν L ν R and g ν L ν R from bearing effects on the dynamics of b → sνν. The current upper bounds on the branching ratios are and supported by the SM results B SM K = BR(B + → K + νν) = (3.98 ± 0.43 ± 0.19) × 10 −6 B SM K * = BR(B + → K * + νν) = (9.19 ± 0.86 ± 0.50) × 10 −6 (25) have been used in Ref. [53] to derive the following 90% C.L. upper bounds We consider the impact of these constraints in the numerical analysis presented in Sec. 4.

The R K ( * ) anomaly
We now turn our attention to the R K ( * ) anomaly, discussing its experimental status and showing how it can be addressed in the present framework.

Experimental status and effective Lagrangian
Another longstanding anomaly highlighted by B physics experiments concerns the neutral current transition b → s + − . More in detail, the LHCb experiment has found anomalous values of and reporting, respectively [8,9], where q 2 is the invariant mass of the final state di-lepton system.
The corresponding SM predictions 3 , suppressed by the GIM mechanism, amount to and giving consequently rise to a discrepancy with a significance of about 5σ [54] depending on the details of the fit.
Global fits of the anomalies presently converge on a preferred sets of Wilson coefficients [20,[54][55][56], with the highest pull given by C 9 − C SM 9 −1.21 or C 9 − C SM 9 = − C 10 − C SM 10 −0.67 that sets an almost 5σ discrepancy with respect to the SM predictions.

Explaining R K ( * ) in the 3HDM
The explanation of the R K ( * ) anomaly in our 3HDM reflects, for the most part, the steps required to account for R D ( * ) one. The possible role of the scalar operators is once again strongly constrained, in this case by observations of B s → µ + µ − . We therefore seek an enhancement of the V − A quark current, along with a suppression of scalar Wilson coefficient, in order to fit the current measurements.
The explanations of the two anomalies share a common origin in the couplings of the quark sector fc L t R and gb L t R . The additional source of lepton flavor universality violation is here implemented through the interactions with the RHNs and muons where Figure 6: The diagram responsible for reproducing the R K ( * ) anomaly within the considered 3HDM.
Our choice of couplings for the charged scalar fields results in vectorial Wilson coefficients that obey C 9 = −C 10 , with the contribution of the diagram in Fig. 6 amounting to where D dd00 is the same 4-points scalar integral given in eq. (18). Notice that the t RsL coupling is obtained from the t RcL one via the CKM element V * cs .

Results
We gather here the results obtained for the R D ( * ) and R K ( * ) anomalies through the analysis detailed in the previous sections, discuss their compatibility within the present framework and remark on the possible impact of new determinations of / on our conclusions.

Numerical analysis
As for the anomalies, the joint effect of the tree-level scalar contribution and of the operator [cγ µ P L b] ¯ γ µ P L ν , generated by the same particles at the one-loop level, allows access to the 2σ and 1σ R D ( * ) regions for values of the Yukawa sector well within the perturbative regime. Remarkably, once supplemented with an extra coupling to muons, we find that the same interactions between quarks and the new scalar doublets allow also to explain the R K ( * ) signal. The result concerning the R D ( * ) anomaly is illustrated in isolation in Fig 7. Here we show four benchmark points in the space of Wilson coefficients obtained through the contributions detailed in Sec. 2. The points in red comply with the 10% lifetime bounds on B c decay, while the blue ones refer to the corresponding 30% limit. The areas shaded in red, orange and yellow indicate the 68%, 95% and 99% confidence interval indicated by current measurements, respectively. The two panels differs by the considered value of the vector Wilson coefficients C τ V L = −C τ AL , induced by the 3HDM at the one-loop level. As shown in the right panel, larger values of this quantity allow to fit the anomaly with a contribution from C τ P L small enough to comply with the B c lifetime bounds. The benchmark points have been obtained by maximizing the one-loop contributions encoded in the vector and pseudovector operators, which depend on the set of couplings rendered in teal in eq. (13). Considering the collider phenomenology analysis presented in Ref. [29], we set the values of the new quark Yukawa couplings to fc L t R = 0.8 and gb L t R = 0.8 or gb L t R = 1. The magnitude of the RHN and lepton couplings f ν L ν R = g ν L ν R and fτ L τ R = gν L τ R , show in Fig 8, is then obtained by setting the Wilson coefficients to the indicated values.
The subdominant tree-level contribution of C τ SL and C τ P L is subsequently obtained through Eq. 14, in compliance with the B c lifetime bounds. The required magnitude of the involved couplings fb L c R and fc L b R , which do not enter the vector contribution, is presented in Fig 9. Here we let the mass of the new scalar particles vary in the selected range and set fν L τ R ≡ fτ L τ R according to the vector operator. Finally, the R K ( * ) anomaly can be accommodated by setting the remaining independent parameters of the model fν L ν R = gν L ν R , resulting in O 9 = −O 10 = −0.67, as required by global fits of the anomalous signal.

Compatibility of the two anomalies
A possible problem for the simultaneous explanations of the charged and neutral current flavour anomalies arises from the definition of R D ( * ) where the denominator contains an average over muons and electrons. It could be consequently thought that the couplings introduced in eq. (12) to explain the R K ( * ) anomaly (in orange), if sizeable, may dilute the yield of the model to R D ( * ) . However, because the contribution of the diagram in the left panel of Fig. 10 is negligible, the analyses of R D ( * ) and R K ( * ) are essentially uncorrelated within the present framework. This is due to the fact that the the couplings employed in the muon sector to explain the R K ( * ) anomaly are much smaller than the ones entering the expression for R D ( * ) , as shown in the right panel of Fig. 10.

Impact of new / determinations
To conclude our analysis we investigate the robustness of our model with respect to another flavour observable that recently received increasing attention [57,58]: / .
This quantity, which quantifies the direct CP violation in K → ππ decays, has been debated in the literature as a further indication of new physics in flavour measurements. The discussion has been revived after the lattice QCD result obtained by the RBC-UKQCD group [59,60] Re( / ) = 1.38(5.15)(4.59) × 10 −4 , was corroborated by a large N c dual QCD estimate [61,62] of The mutual agreement of these independent estimates of the SM contribution results in a 2σ deviation from the current experimental measure [63][64][65] Re( / ) exp = (16.6 ± 2.3) × 10 −4 .
Presently the origin of the anomaly remains unclear as misdeterminations of the SM contribution alone could explain the mentioned deviation [58]. In light of this ambiguous signal, we discuss below how new determinations of the SM contribution into / could constrain, or be rectified, in the present framework.
The contribution of a charged scalar to the effective d → sg Hamiltonian where the chromomagnetic dipole operators (CMOs) are results in a short distance modification of the SM dynamics which, without introducing further sources of CP violation, alleviates the tension with the current experimental results [66,67].
The determination of the hadronic matrix element for K 0 → ππ in dual QCD [68][69][70][71][72], allows an estimate of order of magnitude needed for the Wilson coefficients in eq. (41) to reproduce the / measurement [62] Re / 8g ∼ − 1.85 × 10 5 GeV × Im C − 8g , where C − 8g is the combination We stress that eq. (43) is obtained under the assumption of a SM-like value for the indirect CP violation in K , enforced in the framework at hand by the same bounds targeting the mass splitting and CP violation of the neutral K system [73]. As previously shown in Ref. [29], these constraints do not significantly limit the parameter space selected by the anomalous signals. With the Yukawa texture in eq. (13), the same interactions behind our explanations of the R D ( * ) and R K ( * ) anomalies yield new contributions to the CMOs that source / . In particular, the propagation of H − 2 in Fig. 11 results in , and, in turn, in value of / comparable in magnitude to the current SM estimate. Therefore, should a revised computation of the SM contribution explain the observed discrepancy, precision measurements of this parameter will certainly offer a new way to probe our scenario.
On the contrary, should the / measurement require the presence of new physics, our framework could offer a possible solution through a minimal modification of the proposed Yukawa texture where the new diagonal down-Yukawa coupling fs L t R results, via the H − 1 contribution shown in Fig. 11, in the term Notice the relevance of the Yukawa term fc L t R , which is therefore pivotal for the explanation of all of the anomaly discussed in the present work and offers a clear signature for the investigation of our proposal.
In Fig. 11 we assess the value of fd L t R required to have Re / 8g = 10 −3 , considering the values of fc L t R and gb L t R used in the benchmark points of the previous section. Clearly, it is sufficient to set fd L t R ∼ 10 −1 to enhance the 3HDM contribution into the CMOs and explain the / anomaly on top of the R D ( * ) and R K ( * ) signals.

Discussion and conclusions
We attempted to demystify the implications of the R D ( * ) and R K ( * ) anomalies in the context of physics beyond the standard model. Although these signals seem to indicate new lepton flavour violating physics at very different scales, we demonstrated that this fact does not necessarily imply the existence of exotic new physics. In fact, the full dynamics responsible for all of the analyzed measurements can still be addressed within a more conventional framework, using model building elements already available in the standard model.
Specifically, we have shown that the present R D ( * ) and R K ( * ) anomalies can be simultaneously explained in a 3 Higgs doublets model extended with right-handed neutrinos. The results plotted in Fig. 7 demonstrate that the model predicts values of R D ( * ) within the 1σ region indicated by the experiments, in accordance with the remaining phenomenological constraints and for perturbative values of the involved Yukawa couplings. The result is achieved owing to the interplay between the loop induced vector operators and tree-level scalar operators in eq. (4), which depend on different sets of parameters. We have furthermore shown that in the considered scheme the R D ( * ) and R K ( * ) anomalies arise from independent interactions of the additional scalar doublets, and consequently found a simultaneous explanation for the two measurements.
The robustness of our results was tested against a further flavour observable, / . Once again, we find that within the proposed 3HDM this anomaly can be modeled independently from the remaining flavour physics signals. Should the current deviation in this parameter be confirmed as an effect of new physics, our proposal would then provide a first framework able to simultaneously explain all of the mentioned anomalies.
A crucial aspect of our work is the inclusion of GeV scale right-handed neutrinos, which allowed for the sizable looplevel contributions needed to explain the analyzed signals. It is remarkable that the presence of these particles in Nature is advocated also in connection to other open problems of contemporary physics [74][75][76], and our framework therefore connects the R D ( * ) and R K ( * ) flavour anomalies to the well-known phenomenology of neutrino masses and to the puzzle of the baryon asymmetry of the Universe. The considered right-handed neutrinos are also a primary target of the next generation beam-dump experiment SHIP, at CERN, which has therefore the potential to directly test our scenario on top of complementary collider searches for new scalar particles.
In the light of the proposed scheme, it is natural that the first signals of new physics would be detected in low-energy flavour experiments rather than in high-energy collider searches and precision observables: because the GeV-scale right-handed neutrinos phenomenology is intrinsically a low-energy phenomenon, the first signals of their presence naturally occur in the loop processes behind the flavour observables under discussion, preceding the potential collider signatures of still undiscovered scalar particles.