Linking ${ {\boldsymbol R}_{{\boldsymbol K}^{\boldsymbol{(*)}}}}$anomalies to Hubble tension via a single right-handed neutrino

In rare semi-leptonic B-meson decays, there is a series of long-standing deviations between standard model (SM) predictions and the LHCb measurements [1–6]. In particular, the ratios $ R_{K^{(\ast)}} $, which are defined as

are predicted to be $ R_{K^{(\ast)}}^{\rm{SM}} = 1.00 \pm 0.01 $ in the region $ 1.1 \le q^2 \le 6\, \mathrm{GeV}^2 $, with $ q^2 $ as the dilepton invariant mass squared, within the SM [7–10], whereas the LHCb measurements in both 2017 [3] and 2019 [4] exhibited a deviation at the $ \sim2.5\sigma $ level. Strikingly, the latest update of the LHCb measurement [5], with

has pushed the deviation to reach the level of $ 3.1\sigma $ owing to the reduced experimental uncertainties. This strongly hints at new physics (NP) beyond the SM that violates lepton-flavor universality (LFU).

$ R_{K^{(\ast)}} $ anomalies have prompted many NP proposals under extensive and intensive investigations (see, for example, the recent reviews [11, 12] and references therein). In particular, the two-Higgs-doublet model (2HDM) extended with right-handed neutrinos [13–16] is an interesting NP candidate because it can connect the intriguing LFU violation with neutrino masses – another big mystery in contemporary particle physics. Thus far, either Dirac [13] or Majorana [14–16] neutrinos have been considered to address $ R_{K^{(\ast)}} $ anomalies. In these scenarios, the dominant NP contribution to $ R_{K^{(\ast)}} $ anomalies is assisted by the right-handed neutrinos running in the box diagrams that are most relevant for the $ b\to s \mu^+\mu^- $ transition, and the resulting LFU-violating Wilson coefficients are predicted in the direction $ C_{9\mu}^{\mathrm{NP}}=-C_{10\mu}^{\mathrm{NP}} $, which implies a left-handed NP effect in the muon sector and is persistently favored by updated global fits following Eq. (2) [17–29]. However, as noted in Ref. [15], the LFU-conserving Wilson coefficient $ C_{10\ell, Z}^\mathrm{NP} $ resulting from Z-penguin diagrams could contribute as much as the LFU-violating ones. Even though the LFU-conserving contribution $ C_{10\ell, Z}^\mathrm{NP} $ cannot explain $ R_{K^{(\ast)}} $ anomalies alone, its comparable contribution can affect how $ R_{K^{(\ast)}} $ anomalies are numerically addressed in the direction $ C_{9\mu}^{\mathrm{NP}}=-C_{10\mu}^{\mathrm{NP}} $. Therefore, theoretical NP models with comparable $ C_{9\mu}^{\mathrm{NP}}=-C_{10\mu}^{\mathrm{NP}} $ and $ C_{10\ell, Z}^\mathrm{NP} $ should match a two-parameter global fit following Eq. (2) [24], which have unfortunately been neglected in Refs. [13, 14, 30].

Many previous studies have focused on heavy Majorana neutrinos for the $ R_{K^{(\ast)}} $ resolution. Nevertheless, as found in Ref. [14], the solution of $ R_{K^{(\ast)}} $ anomalies in the direction $ C_{9\mu}^{\mathrm{NP}}=-C_{10\mu}^{\mathrm{NP}} $ is insensitive to Majorana neutrino masses below the electroweak scale. This implies that the difference between heavy and light Majorana neutrinos cannot be simply distinguished by the $ R_{K^{(\ast)}} $ resolution. Besides Majorana neutrinos, topics with Dirac neutrinos have also received increasing attention in recent years, especially in connection with the phenomenologies [31–36] of big-bang nucleosynthesis (BBN) and cosmic microwave background (CMB), as well as the baryon asymmetry of the Universe [37, 38]. Dirac neutrino effects in the $ b\to s \mu^+\mu^- $ process have been found, as noted in Ref. [13]. However, it was found that because $ \mathcal{O}(1) $ Dirac neutrino Yukawa couplings are generically required to explain $ R_{K^{(\ast)}} $ anomalies, thermalized right-handed Dirac neutrinos with such large couplings in the early Universe would cause an undesired shift in the effective neutrino number, $ \Delta N_{\rm eff}=N_{\rm eff}-N_{\rm eff}^{\rm SM} $, at the BBN and CMB epochs, where the SM prediction reads $ N_{\rm eff}^{\rm SM} = $ 3.044– 3.045 [39–45]. Nevertheless, the above conclusion depends crucially on how many right-handed Dirac neutrinos are thermalized in the early Universe and the decoupling temperature of thermalized neutrinos, which can result in different levels of the $ \Delta N_{\rm eff} $ shift.

In addition to the extensive investigations of the heavy nature of Majorana neutrinos used to address $ R_{K^{(\ast)}} $ anomalies, it is also interesting to consider situations in which eV-scale Majorana neutrinos are included. In this paper, we consider NP effects on $ R_{K^{(\ast)}} $ anomalies from either an eV-scale right-handed Majorana or a right-handed Dirac neutrino. Although the NP effects arising from these two cases are indistinguishable in terms of the $ R_{K^{(\ast)}} $ resolution alone, their impacts on the early Universe are different in generating an observable $ \Delta N_{\rm eff} $ shift owing to the spinor nature of the neutrinos involved, that is, the Majorana spinor for the former and the Weyl spinor for the latter. Therefore, it becomes possible to distinguish these two solutions via the observation of different $ \Delta N_{\rm eff} $ shifts in the cosmic regime.

Noticeably, the extra radiation that generates a significant $ \Delta N_{\rm eff} $ shift is one of the simplest candidates to mitigate Hubble ($ H_0 $) tension (see, for example, Refs. [46–50] for the latest reviews), which signifies a notorious discrepancy between the local measurements of the present-day Hubble parameter from the SH0ES collaboration [51–53] (based on the publication years, the three references are dubbed R18, R19, and R20) and the Planck CMB inferred value under the standard ΛCDM baseline [54],

The $ H_0 $ tension is further worsened by the updated SH0ES measurements (R21) [55], with

enhancing the deviation from the Planck 2018 data to $ 5\sigma $. A plethora of investigations have introduced a shift in the effective neutrino number, $ \Delta N_{\rm eff}\simeq 1.0 $, to address $ H_0 $ tension [56–59]. As illustrated in Ref. [58], an extra free $ N_{\rm eff} $ beyond the original six ΛCDM parameters can cause a genuine shift in the central value of $ H_0 $ from Planck measurements, and $ H_0 $ tension can be relieved with $ N_{\rm eff} \approx 3.95 $. Here, $ N_{\rm eff} $ serves as an NP source to shift the ΛCDM predictions inferred from CMB, BAO, and Pantheon Supernovae Type-Ia data so that they are in agreement with the local $ H_0 $ measurements. As opposed to estimating $ N_{\rm eff} $ simply by combining the high-redshift measurements with local $ H_0 $ data (like in many other studies), the data-analyzing method proposed in Ref. [58] opens a new avenue to ease $ H_0 $ tension. However, such a large shift is disfavored by high-$ \ell $ Planck CMB polarization measurements [54, 60–62]. More recent analyses have instead showed that a shift of $ 0.2<\Delta N_{\rm eff}<0.6 $ is able to ease $ H_0 $ tension. For instance, Ref. [63] points out two possible regimes with/without BBN data,

in which the SH0ES 2019 measurements (R19) [52] are included. These patterns are also consistent with that observed in Ref. [64], where an additional electron-type lepton asymmetry $ \xi_e $ in the neutrino sector is introduced, giving

with a larger central value of $ N_{\rm eff} $ than from Eq. (6). Intriguingly, the introduction of lepton asymmetry is also supported by the very recently probed anomaly in the helium-4 abundance [65], which results in [65–67]

Note that all the central values of $ N_{\rm eff} $ obtained above are larger than the previous CMB+BBN result [68], $ N_{\rm eff}=2.843\pm 0.154 $. Therefore, it can be inferred from Eqs. (5)–(8) that an increased $ N_{\rm eff} $ will be helpful to mitigate $ H_0 $ tension, though an updated analysis of the combined dataset from CMB+BAO+Pantheon+R21+BBN is currently not available. Other possible patterns, such as the extra radiation in the presence of additional non-free-streaming degrees of freedom (d.o.f) [69, 70], also found that comparable $ N_{\rm eff} $ values are favored to ease $ H_0 $ tension. Recently, the mitigation of $ H_0 $ tension with $ \Delta N_{\rm eff}\simeq \mathcal{O}(0.5) $ has been studied in several explicit models [62, 71–73].

The results above suggest that a full resolution of $ H_0 $ tension could be a result of multidisciplinary interplay, in which the extra radiation serves a fractional but important role. In relating the observed anomalies in the particle physics domain, it is compelling to consider the situation where the underlying mechanism for the $ \Delta N_{\rm eff} $ shift is naturally provided by the $ R_{K^{(\ast)}} $ resolution via an eV-scale right-handed Majorana or a right-handed Dirac neutrino, which motivates our present study. In this paper, we show that such a connection can indeed be realized in a flavor-specific 2HDM framework, where only one right-handed Majorana or Dirac neutrino has significant interactions with the extra Higgs bosons present in the model.

The paper is organized as follows. We begin in Sec. II with a description of the framework, dubbed $ t\nu $2HDM, and then consider in Sec. III the most relevant constraints from low-energy flavor physics, the perturbative unitarity condition, and LHC direct searches. In Sec. IV, we discuss the NP contributions to $ R_{K^{(\ast)}} $ anomalies and the mitigation of $ H_0 $ tension. Then, we present in Sec. V our detailed numerical analyses of the viable parameter space for the $ R_{K^{(\ast)}} $ resolution, as well as the correlation between $ R_{K^{(\ast)}} $ anomalies and $ H_0 $ tension. Conclusions are presented in Sec. VI.

The 2HDM is a simple extension of the SM achieved by adding a second Higgs doublet to the SM particle content [74, 75]. Any specific 2HDM framework is characterized by its Yukawa interactions and scalar potential, both of which can be specified by either symmetry backgrounds or purely phenomenological considerations. For our purpose in addressing $ R_{K^{(\ast)}} $ anomalies with a link to $ H_0 $ tension, we follow a data-driven approach.

The two Higgs doublets $ H_{1, 2} $ in the model are constructed in the so-called Higgs basis [75, 76] as

where the vacuum expectation value $v\simeq 246~ \mathrm{GeV}$ is responsible for generating the fermion and gauge-boson masses, and $ G^{+, 0} $ are the Goldstone bosons. Here, we assume a CP-conserving Higgs potential [75]. Then, $ H^+ $ and A are the physical charged and neutral pseudo-scalar Higgs bosons, respectively, whereas the neutral scalars $ \phi_{1, 2} $ are the superposition of the two mass eigenstates $ H^0 $ and h, which can generally be written as

with the mixing angle θ determined completely by the parameters in the Higgs potential. Given that $ \cos\theta\approx 1 $ is favored by current LHC data on various SM-like Higgs signals (see, for example, Refs. [77, 78] for recent global fits of the 2HDMs), which corresponds to the so-called alignment limit, we consider the case where $ H_1 $ is the SM Higgs doublet such that h corresponds to the observed Higgs boson [79, 80], whereas $ H_2 $ is the NP doublet with $ H^0 $ corresponding to the extra physical neutral scalar.

In the alignment limit, the Higgs potential $ V(H_1, H_2) $ can be readily constructed in terms of the free potential parameters governing the Higgs mass spectrum. In principle, these free parameters receive various theoretical and phenomenological constraints, such as vacuum stability, perturbative unitarity, electroweak precision tests, and collider direct detection [77, 78]. Nevertheless, the mass spectrum of the physical states $ H^+ $, $ H^0 $, and A is still undetermined by current LHC direct searches. In particular, a quasi-degenerate Higgs mass spectrum,

still remains a possible regime and is covered here.

To address $ R_{K^{(\ast)}} $ anomalies with a link to $ H_0 $ tension, we consider the following Yukawa interactions:

where $ \mathcal{L}_Y^{\rm{SM}}(H_1) $ denotes the SM Yukawa Lagrangian associated with the Higgs doublet $ H_1 $, whereas $ \mathcal{L}_Y (H_2) $ encodes the NP interactions related to the second Higgs doublet $ H_2 $, with $\tilde{H}_2= {\rm i}\sigma_2 H_2^\ast$, and $ \sigma_2 $ as the second Pauli matrix. The left-handed fermion doublets $ \bar Q_L $ and $ \bar E_L $ are specified as

respectively, where all the chiral fermions $ f_{L, R} $ ($ f=u, d, e, \nu $) are physical fields, with V and U representing the Cabibbo-Kobayashi-Maskawa (CKM) and Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrices, respectively. For the Yukawa matrices $ X_{u, \nu} $, we propose the following phenomenologically viable structure:

where $ \kappa_{t, \nu} $ are the only nonzero real effective couplings, and the flavor index s characterizes the one-flavor right-handed neutrino that couples to the muon lepton in the charged scalar current. It should be noted that the explicit right-handed neutrino flavor is irrelevant here and will simply be denoted as $ \tilde{\nu}_{R} $ hereafter.

Our proposal of Eq. (15) comes from various data-driven considerations. In the quark sector, Eq. (13) together with Eq. (15) would induce only neutral scalar currents associated with the top quark, and the charged scalar interactions,

have only significant effects on the third generation of quarks owing to the hierarchical structure of the CKM matrix. These patterns comply with the current observation that only significant NP contributions are allowed in the third generation and the flavor-changing neutral scalar currents are severely constrained by experimental data [81–83]. In the lepton sector, however, Eqs. (13) and (15) would indicate that there are only neutral scalar currents in the neutrino sector, whereas the charged scalar interaction is only stimulated by the one-flavor right-handed neutrino $ \tilde{\nu}_{R} $ that has a significant coupling to the muon lepton, namely,

Such particular patterns closely follow the tight bounds from the charged lepton-flavor violating processes $ \ell_i\to \ell_j \gamma $ mediated by the right-handed neutrino at the loop level [14] and the muon decay $ \mu\to e \nu \bar{\nu} $ mediated by the charged Higgs at the tree level. Furthermore, the reason for allowing only one rather than two or three flavors of right-handed neutrinos to interact with the muon lepton comes from the observation that if more than one right-handed neutrino has significant couplings to account for $ R_{K^{(\ast)}} $ anomalies, the resulting parameter space will readily force them to establish thermal equilibrium in the early Universe and thus generate an unacceptably large $ \Delta N_{\rm eff} $ shift [13, 32], as confirmed later in this paper. Finally, it should be emphasized that, because we are interested in the connection between $ R_{K^{(\ast)}} $ anomalies and $ H_0 $ tension via a minimal setup, other couplings not considered in $ X_{u, \nu} $ do not necessarily need to vanish, but rather they signify the meaning of phenomenological smallness in their own right. Furthermore, we are not concerned with the symmetry underlying such a flavor-specific Yukawa structure given by Eq. (15), though interesting possibilities, such as Branco-Grimus-Lavoura-based scenarios [84] and mass-powered-like textures [85], may deserve further exploitation.

The above considerations result in our flavor-specific 2HDM framework that points toward significant NP effects associated with the top quark t and one-flavor right-handed neutrino $ \tilde{\nu}_R $ and will therefore be dubbed the $ t\nu $2HDM hereafter. As mentioned in Sec. I, the neutrino nature, being of the Majorana or Dirac type, is unspecified by the $ R_{K^{(\ast)}} $ resolution alone. However, when $ \tilde{\nu}_R $ is an eV-scale Majorana neutrino, its impact on the $ \Delta N_{\rm eff} $ shift will be different from that with a Dirac neutrino, especially when the shift is linked to the mitigation of $ H_0 $ tension. Furthermore, if $ \tilde{\nu}_R $ is of the Majorana type, it can be embedded into the seesaw mechanism (see, for example, Refs. [86, 87] for recent comprehensive reviews), where two more right-handed Majorana neutrinos are introduced with the presence of a Majorana mass term,

where $ (\nu_R)^c=C\overline {\nu_R}^T $, and $C={\rm i} \gamma^2\gamma^0$ is the charge conjugation matrix. Then, the active neutrino masses are generated via the seesaw formula,

in the basis where the symmetric matrix $ M_R $ is already diagonal. Here, $ Y_\nu $ is the neutrino Yukawa matrix from $ \mathcal{L}_Y^{\rm{SM}}(H_1) $ in Eq. (12), that is,

encoding active-sterile neutrino mass mixing and mixing-induced interactions via [88]

An important observation arises if one of the sterile neutrino eigenstates is at the eV-scale, for example, $ M_{R, 1}\simeq \mathcal{O}(\mathrm{eV}) $. To avoid the constraints on active-sterile neutrino mixing [89], particularly for an eV-scale sterile neutrino [90], the first column of $ Y_{\nu} $ must be strongly suppressed. Otherwise, $ W_{\nu, i1} $ would be enhanced by a factor of $ v/M_{R, 1}\simeq 10^{11} $. In the asymptotically safe limit, $ Y_{\nu, i1}=0 $, the active neutrino mass matrix from Eq. (19) would be of rank two, making the lightest active neutrino massless. Therefore, if the eV-scale $ \tilde{\nu}_R $ belongs to the lightest sterile neutrino in the seesaw mechanism, the lightest active neutrino in the $ 3\nu $ oscillation paradigm [91] is essentially massless. Further constraints on $ W_{\nu, ij} $ are not discussed in the following because the NP effects under consideration in this paper do not rely on the neutrino Yukawa matrix $ Y_{\nu} $. If $ \tilde \nu_R $ is of the Dirac type, however, the Dirac neutrino mass may also be generated via Eq. (20) but with the absence of the Majorana mass term given by Eq. (18). In either case, Eq. (13) will encode all the NP interactions concerned in this paper.

In the following sections, we consider important constraints from flavor and collider physics, as well as perturbative unitarity on the $ t\nu $2HDM framework, which is only characterized by the three free parameters $\kappa_t, ~ \kappa_\nu, ~ m_S, $ and show the viable parameter space in addressing $ R_{K^{(\ast)}} $ anomalies. We further show that the resulting favored parameter space induces the shift $ \Delta N_{\rm eff}\simeq 1.0 $ in the Dirac neutrino case and $ \Delta N_{\rm eff}\simeq 0.5 $ in the Majorana case.

The $ t\nu $2HDM framework indicates significant NP effects associated with the third generation of quarks and the muon lepton. In this section, we discuss the most relevant constraints on the model from low-energy flavor physics, the perturbative unitarity condition, and LHC direct searches.

In the framework of low-energy effective field theory, the effective Hamiltonian governing the radiative $ b\to s\gamma $ decay at the scale $ \mu_b \simeq \mathcal{O}(m_b) $ reads as

where $G_{\rm F}=1/(\sqrt{2}v^2)$ is the Fermi constant, and the terms proportional to $ V_{us}^\ast V_{ub} $ are neglected in view of $ |V_{us}^\ast V_{ub}/V_{ts}^\ast V_{tb}|<0.02 $. The explicit expressions of the current-current ($ \mathcal{O}_{1, 2} $) and QCD-penguin ($ \mathcal{O}_{3-6} $) operators can be found in, for example, Refs. [92–95], whereas the magnetic dipole operators are defined, respectively, by

where $ P_{R, L}=(1\pm\gamma_5)/2 $ are the right- and left-handed chirality projectors.

In the $ t\nu $2HDM framework, the NP contributions to the Wilson coefficients $ C_{1-6} $ are absent, and their contributions to the primed dipole coefficients $ C_{7\gamma, 8g}^\prime $ are suppressed by the ratio $ m_s/m_b $. As a consequence, the dominant NP influence on the $ b\to s\gamma $ transition stems from the unprimed $ C_{7\gamma} $ and $ C_{8g} $. After a direct evaluation of the one-loop penguin diagrams with the charged Higgs running in the loop, as shown in Fig. 1, we can obtain the NP Wilson coefficients at the matching scale $ \mu_S\simeq \mathcal{O}(m_S) $ [93–95],

Figure 1.  Relevant NP photon- (the first two) and gluon-penguin (the last) diagrams contributing to the inclusive radiative $ \bar{B}\to X_s\gamma $ decay.

where the scalar functions are defined, respectively, by

with $ z_t \equiv m_t^2/m_S^2 $, and $ m_t $ is the top-quark $ \overline {\rm{MS}} $ mass.

To evaluate the NP contributions to the branching ratio $ \mathcal{B}(\bar{B}\to X_s\gamma) $, we must run the Wilson coefficients from the matching scale $ \mu_S $ down to the low-energy scale $ \mu_b $ [96, 97]. Generically, the Wilson coefficient $ C_{7\gamma}(\mu_b) $ can be divided into two parts,

which contributes to $ \mathcal{B}(\bar{B}\to X_s\gamma) $ with a photon-energy cutoff $E_\gamma < 1.6\,\ \mathrm{GeV}$ via [96, 97]

Here, the overall factor reads numerically as $ R = 2.47 \times 10^{-3} $, and we neglect the small non-perturbative effect. The SM contribution $ C_{7\gamma}^{\rm{SM}}(\mu_b) $ has been calculated up to the next-to-next-to-leading order in QCD [98–100], and the resulting numerical value reads as [101, 102]

whereas the NP part $ C_{7\gamma}^{\rm{NP}}(\mu_b) $ is given by

where $ C_{7\gamma}^{\rm{NP}}(\mu_S) $ and $ C_{8g}^{\rm{NP}} (\mu_S) $ are already given by Eq. (24), and the magic numbers are evaluated to be $ \kappa_7 =0.457 $ and $ \kappa_8 = 0.125 $ at the NP scale $ \mu_S \sim \mathcal{O}(1)\, \mathrm{TeV} $.

By comparing the theoretical prediction given by Eq. (28) with the current world-averaged experimental data [103],

we can set bounds on the NP Wilson coefficients presented in Eq. (24), and the allowed parameter space $ (\kappa_t, m_S) $ can therefore be extracted. In Sec. V, we apply the $ \mathcal{B}(\bar{B}\to X_s\gamma) $ constraint within a $ 1\sigma $ uncertainty.

Next, we turn our attention to the mass differences $ \Delta M_{d, s} $ describing the strength of $ B^0_{d, s}-\bar{B}^0_{d, s} $ mixings. The theoretical description of $ B^0_{d, s}-\bar{B}^0_{d, s} $ mixings can be realized in terms of the low-energy effective Hamiltonian

where $ m_W $ is the W-boson mass, and $ q=d(s) $ for the neutral $ B_{d(s)} $ meson. The explicit expressions of the four-quark operators can be found in, for example, Refs. [104, 105].

In both the SM and $ t\nu $2HDM framework, the only significant Wilson coefficient responsible for neutral B-meson mixing originates from $ C_{1q} (\mu_b) $, which corresponds to the four-quark operator

where the Greek letters α and β denote the quark color indices. The mass difference of neutral B-meson mixing can be expressed in terms of the off-diagonal matrix element, $ \Delta M_q = 2|\mathcal{M}_{12}^q| $, with the latter given by [104, 105]

Here, the hadronic matrix element $ \langle \bar{B}_q^0| \mathcal{Q}_{1q} |B_q^0 \rangle $ encodes the non-perturbative QCD effect, whereas the perturbative contribution is absorbed into the short-distance Wilson coefficient $ C_{1q} (\mu_b) $. Normalizing NP to the SM contribution, we can parameterize the theoretical prediction of $ \Delta M_q $ as

where $ \Delta M_q^{\rm{SM}} $ denotes the SM contribution. For the NP scale of $ \mathcal{O}(1)\, \mathrm{TeV} $ considered in this paper, the correction $ \Delta^{\rm{NP}}_q $ is given by

where $ U^{(0)}(\mu_i, \mu_j) $ represents the leading-order QCD evolution function from the high-scale $ \mu_j $ to low-scale $ \mu_i $ [104], and $ C_{1q}^{\rm{SM}}(\mu_W) $ is the SM Wilson coefficient evaluated at $ \mu_W \simeq \mathcal{O}(m_W) $. Here, we take into account the threshold effect when evolving across the top-quark mass scale $ \mu_t\simeq \mathcal{O}(m_t) $ [104]. The NP contribution $ C_{1q}^{\rm{NP}}(\mu_S) $ is obtained by evaluating the box diagrams shown in Fig. 2, leading to

Figure 2.  Relevant NP box diagrams responsible for $ B^0_{q}-\bar{B}^0_{q} $ mixing in the 't Hooft-Feynman gauge.

in which the different parts are given, respectively, as

where $ z_W \equiv m_W^2/m_S^2 $, and we introduce the following scalar functions:

The current world-averaged experimental measurements are [106]

both of which carry considerably smaller uncertainties than those of the corresponding theoretical predictions [107–110]. In particular, based on the bag parameters calculated in full four-flavor lattice QCD for the first time, the HPQCD collaboration found that [108]

in which the central value of $ \Delta M_d^{\rm{SM}} $ is larger than the experimental data. This in turn implies a discrepancy for the ratio $ \Delta M_d/\Delta M_s $ at $ \sim 1.7\sigma $. However, an earlier computation based on the most accurate numerical inputs at that time found that [107]

the central value of which is $ \sim 1.8\sigma $ above the experimental one given by Eq. (44). Such a difference has profound implications for NP models that predict sizable positive contributions to $ B^0_{s}-\bar{B}^0_{s} $ mixing [107]. While the discrepancies observed in $ \Delta M_{d, s} $ are not yet conclusive owing to the large theoretical uncertainties, it is interesting to note that an excess over the SM predictions cannot be reconciled with the $ R_{K^{(\ast)}} $ resolution in the $ t\nu $2HDM framework because NP effects on $ \Delta M_{d, s} $ are always positive, as shown from Eq. (37). Therefore, if confirmed with more precise experimental measurements and theoretical predictions, the discrepancy will entail additional NP sources beyond the minimal setup considered in this study.

In view of the above observations, we apply in this study the HPQCD results for $ \Delta M_{d, s} $ given by Eq. (45) as constraints but vary the uncertainties within $ 3\sigma $ conservatively. We would like to emphasize again that the constraining power from $ \Delta M_{d, s} $ can be highly efficient only when the theoretical uncertainties from, for example, B-meson decay constants, bag parameters, and CKM elements are significantly reduced [107].

Besides the B-meson observables discussed above, the $ t\nu $2HDM also has an impact on K-meson observables, such as the branching ratios of the $ K_{S, L}\to\mu^+\mu^- $ decays as well as the mass difference $ \Delta M_K $ and the $ \epsilon_K $ parameter of $ K^0 $-$ \bar{K}^0 $ mixing. However, because kaons are composed of two light quarks, that is, the up (down) and strange quarks, whereas NP interactions in the quark sector within our framework always connect with the top quark (see Eq. (15)), their leading contributions to K decays and mixing must first stem from one-loop diagrams with the top quark and charged Higgs running in the loop. This implies that the impact of NP on K-meson observables is suppressed by both the loop factor and these heavy particle masses, as well as the CKM entries involved.

We explicitly evaluate the short-distance NP contributions to the branching ratios of the $ K_{L, S}\to \mu^+\mu^- $ decays and find that they only result in a negligible effect on the branching ratios of the $ K_{L, S}\to \mu^+\mu^- $ decays, especially when the sign of the long-distance contribution is chosen to be destructive with the short-distance part [111–113]. We also check if the resulting parameter space of the $ t\nu $2HDM complies with the constraint from $ K^0-\bar{K}^0 $ mixing. To this end, fixing the free parameters at a typical benchmark point $ (\kappa_t, m_S) \sim (1.0, 1000\, {\rm GeV}) $, we numerically find a significantly weaker impact on $ K^0-\bar{K}^0 $ mixing compared with that obtained through a global fit study [81]. Thus, from these observations, we may safely conclude that K-meson observables do not place any significant constraints on the $ t\nu $2HDM compared with those obtained from B-meson observables.

As a consequence, we do not show the constraints from K-meson observables in the numerical analysis below.

Let us now consider the constraints from the LFU ratios of the di-lepton decays of Z and W gauge bosons, $ \Gamma(Z\to \mu^+\mu^-)/\Gamma(Z\to \ell^+\ell^-) $ and $ \Gamma(W\to \mu \bar{\nu})/\Gamma(W\to \ell \bar{\nu}) $, where $ \ell=e $ or τ. For both of these two cases, by encoding one-loop NP corrections into the renormalized effective vertex in the on-shell scheme, we can readily derive the NP contributions to these LFU ratios.

For Z-boson decays, the LFU ratio $ R_{\mu \ell}^Z $ can be parameterized as

in the vanishing lepton mass limit. Here, $ R_{\mu \ell}^{Z, \rm{SM}} $ is the SM contribution, and the SM couplings are given by $ g_{V, Z}^{\rm{SM}} = -1/2 + 2 s_W^2 $ and $ g_{A, Z}^{\rm{SM}} = -1/2 $, with $ s^2_W\equiv \sin^2 \theta_W \simeq 0.23 $ [114]. It should be noted that the NP contribution to the electron/tauon mode is absent in view of the flavor-specific Yukawa structure given by Eq. (15). Given that the charged Higgs only couples to the left-handed muon (cf. Eq. (17)), whereas neutral scalars do not interact with the muon (cf. Eq. (13)), the NP contribution to the LFU ratio originates solely from the $ H^+ $-mediated loop diagram, as shown by the left Feynman diagram in Fig. 3. We explicitly arrive at

Figure 3.  NP contributions to the di-lepton decays $ Z\to \mu^+\mu^- $ (left) and $ W\to \mu \bar{\nu} $ (right).

where $ m_Z $ is the Z-boson mass, and $ c_{2W} = 1-2s_W^2 $. Additionally, in the quasi-degenerate limit for the Higgs mass spectrum, as given by Eq. (11), the NP contributions to the decay $ Z\to \mu^+\mu^- $ from the two neutral Higgs bosons H and A cancel out to a large extent, leaving the dominant NP effect from the left Feynman diagram shown in Fig. 3.

For the W-boson decays, the NP effect arises from the right Feynman diagram shown in Fig. 3. The resulting expression for the LFU ratio $ R_{\mu \ell}^W $ in the vanishing lepton mass limit can be analogously obtained by replacing Z with W in Eq. (47). The corresponding effective couplings are now given by $ g_{V, W}^{\rm{SM}} = g_{A, W}^{\rm{SM}} =1/2 $, and

Note that in deriving the above equation, we make use of the quasi-degenerate Higgs mass spectrum given by Eq. (11).

Then, by comparing the theoretical predictions with the experimental data [114, 115]

we can extract the bounds on the NP parameter space $ (\kappa_\nu, m_S) $. To this end, we take the constraints as the experimental data within $ 1\sigma $ uncertainties.

In addition to the severe constraints from low-energy flavor physics, theoretical considerations, such as the bounded-from-below limit on the scalar potential and the perturbative unitarity condition of high-energy scattering amplitudes (see, for example, Ref. [75] and references therein), may also render tight bounds on NP parameter space. Here, we consider the vital requirement of perturbative unitarity for the Yukawa sector [116].

Generally, the perturbative unitarity bounds can be derived using the so-called partial wave expansion approach [117]. For the case of $ 2\to 2 $ scattering processes in the high-energy massless limit, the partial waves $ a_{fi}^J $ with fixed total angular momentum J are explicitly defined by [117]

where $ d_{\mu_i\mu_f}^J (\theta) $ are small Wigner d-functions, with $ \mu_i=\lambda_{i_1}-\lambda_{i_2} $ and $ \mu_f=\lambda_{f_1}-\lambda_{f_2} $ characterizing the total helicities of the initial and final states, respectively, and $ \mathcal{T}_{fi} (\sqrt{s}, \cos \theta) $ are the invariant scattering amplitudes $ (2\pi)^4\delta^{(4)}((p_{i_1}+p_{i_2})-(p_{f_1}+p_{f_2}))i\mathcal{T}_{fi} (\sqrt{s}, \cos \theta)=\langle f|S-1|i \rangle $, which are related to the S matrix via $S=1+ {i}\mathcal{T}$. Here, θ is the polar scattering angle in the center-of-mass frame, and $ \sqrt{s} $ is the center-of-mass energy. By focusing on the elastic channels, that is, $ |i \rangle =|f \rangle $, corresponding to forward scattering with the same spin and internal variables in the initial- and final-state configurations, and restricting the sum over the intermediate states only to two-particle states, we can obtain from the unitarity condition of the S matrix, $ S^\dagger S=1 $, the following reliable bounds on the tree-level partial-wave scattering matrices [116]:

which give a reasonable indication of the range of validity of the perturbative expansion.

To extract the best perturbative unitarity bounds from Eq. (53), we must then identify the optimal elastic channels. To this end, we must first obtain concrete expressions for the scattering amplitudes $ \mathcal{T}_{fi} (\sqrt{s}, \cos \theta) $, which depend on the definite Yukawa structure and scalar potential as well as the underlying symmetry properties of the model under consideration. Using $ \mathcal{T}_{fi} (\sqrt{s}, \cos \theta) $, it is then straightforward to obtain the perturbative unitarity bound for each entry of Eq. (53) by performing the convolution with the Wigner d-functions and integration over the polar angle θ (cf. Eq.(52)) and then finding the eigenvalues of the partial-wave matrices $ a_{fi}^J $. For generic fermionic Yukawa interactions, owing to the presence of different spin states in the scattering processes, this can be most efficiently achieved in the Jacob-Wick formalism [117]. However, the traditional method for calculating $ \mathcal{T}_{fi} (\sqrt{s}, \cos \theta) $ relies on computing all the matrix entries, which becomes very involved and highly inefficient when the transition matrix has large dimensions. Recently, it was noticed that the determination of perturbative unitarity bounds in this case can be simplified by decomposing each scattering amplitude with different J into a Lorentz part that depends only on the spin and helicity of the fields involved and a group-theoretical part that depends only on their symmetry quantum numbers [116]. The only complication in the method is then attributed to the calculation of the symmetry factors, while the Lorentz parts are universal for different group structures [116].

To obtain the perturbative unitarity bounds on the Yukawa parameters $ \kappa_\nu $ and $ \kappa_t $ of the $ t\nu $2HDM, we employ the results derived in Ref. [116]. For the lepton part, which is characterized by the SM gauge group $S U(2)_L \times U(1)_Y$, the most stringent bound originates from the P-wave amplitude, that is, $ J=1 $, and imposes an upper limit on the muon-related coupling $ \kappa_\nu $.

For the quark part, however, the constraints are different because quarks carry an additional color quantum number under the gauge group $S U(3)_C$. As a consequence, the tightest constraint on the top-related coupling $ \kappa_t $ stems from the S-wave amplitude with $ J=0 $, which leads to

and hence a more stringent bound than on $ \kappa_\nu $.

In the $ t\nu $2HDM framework, as the quasi-degenerate mass regime in the alignment limit is considered, we can see that the $ H^0 $ and A decay modes $ H^0/A\to AZ/ H^0 Z $, $ H^0/A\to H^\pm W^\mp $, and $ H^0\to AA, H^+ H^- $ are all forbidden, and the tree-level triple couplings $ H^0/A-V-V $ (where V denotes one of the gauge vector bosons $ W/Z/\gamma $ and gluons) and $ H^0/A-Z(h)-h $ are absent. This implies that for the heavy scalars concerned in this paper, their dominant decay modes are the top-quark and neutrino pair, whereas the di-boson modes are suppressed by the loop factor and, more importantly, by the mass ratio $ m_t/m_S $ [118]. Therefore, the decay width of the neutral scalars is approximately given by

where $ n_S=3/2 $ for $ S=H^{0} $, and $ n_S=1/2 $ for $ S=A $.

Currently, LHC direct searches for neutral scalar productions have been performed by ATLAS with a center-of-mass energy $ \sqrt{s}=8\, \mathrm{TeV} $ [119] and the CMS collaboration with $ \sqrt{s}=13\, \mathrm{TeV} $ [120] in the channel $ pp\to S\to t\bar t $. In particular, the CMS results set model-independent constraints on the coupling modifiers $ g_{S\bar t t} $ between the scalar S and the top quark:

The exclusion limits on $ g_{S\bar t t} $ can then be translated into the allowed regions of the $ t\nu $2HDM free parameters $ \kappa_{t, \nu} $ and $ m_S $. To this end, we must notice that the exclusion limits set by the CMS collaboration are obtained by assuming a fixed decay width $ \Gamma_S $, with $ \Gamma_S/m_S=[0.5, 25]\% $. However, as inferred from Eq. (56), for $ \kappa_{t, \nu}\sim \mathcal{O}(1) $ and $ m_S\gtrsim 500\, \mathrm{GeV} $, a ratio of $\Gamma_S/m_S\gtrsim (4\%-5\%)$ is obtained. As a consequence, we only apply the two benchmark points $ \Gamma_S/m_S=10\% $ and $ \Gamma_S/m_S=25\% $, selected in Ref. [120], to obtain a rough constraint on $ \kappa_{t, \nu} $ for a fixed scalar mass.

More significant constraints on the model parameters have arisen from LHC direct searches for the charged Higgs performed over the past few years. Both the ATLAS and CMS collaborations have covered several decay channels of the charged Higgs, which are dominated by the $ \tau\nu $ [121, 122] and $ tb $ [123, 124] final states. Recently, it was noticed in Ref. [125] that the $ \mu \nu $ final state can also be an excellent complementary discovery channel of the charged Higgs. However, a comprehensive search for such a channel at the LHC is not yet available, and thus there is no significant bound on NP parameter space from the decay. Specific to the $ t\nu $2HDM framework, the decay modes of the charged Higgs are dominated by the $ tb $ and $ \mu\nu $ final states, and the $ \tau\nu $ mode is suppressed under the flavor-specific Yukawa structure of Eq. (15), with the total decay width given approximately by

where we neglect the bottom-quark and muon masses.

To obtain the viable parameter space of $ (\kappa_{t}, \kappa_\nu, m_S) $, we apply the latest results from ATLAS [123] and CMS [124] with $ \sqrt{s}=13\, \mathrm{TeV} $, where the model-independent exclusion limits on the $ tb $-associated production cross section $ \sigma(pp\to H^\pm t b) $ times the branching fraction $ \mathcal{B}(H^\pm \to tb) $ are obtained for the charged-Higgs mass at $ [0.2,~ 2]\, \mathrm{TeV} $ and $ [0.2,~ 3]\, \mathrm{TeV} $, respectively, although the constraints from the CMS results are weaker than those from the ATLAS searches. For the theoretical prediction in the $ t\nu $2HDM framework, we calculate the cross section $ \sigma(pp\to H^\pm t b) $ using the computer program $ \rm{MadGraph5\_aMC@NLO}$[126], with the charged-Higgs decay width given by Eq. (58).

The right-handed neutrino $ \tilde{\nu}_R $ with its interaction specified by Eq. (17) contributes to the $ b\to s\ell^+\ell^- $ process mainly via the box diagram shown in Fig. 4(a). Its contribution can be described by the effective weak Hamiltonian

Figure 4.  (a): NP box diagram contributing to the $ b \to s \mu^+ \mu^- $ transition, where only one-flavor right-handed neutrino $ \tilde{\nu}_{R} $ participates non-negligibly in the loop. (b): $ \gamma/Z $-mediated NP penguin diagrams contributing to the $ b\to s\ell^+\ell^- $ transition.

where $ \alpha_e=e^2/(4\pi) $ is the fine-structure constant, and the two effective operators are defined, respectively, as

with the corresponding LFU-violating Wilson coefficients $ C_{9\mu}^{\mathrm{NP}} $ and $ C_{10\mu}^{\mathrm{NP}} $ given by

Note that we neglect the neutrino mass in the above formula, and our result is consistent with that obtained in Ref. [14] within the vanishing neutrino mass limit.

In addition to the LFU-violating contribution given by Eq. (61), the flavor-specific Yukawa texture characterized by Eq. (15) also renders a considerable LFU-conserving effect on the $ b\to s \ell^+\ell^- $ transition via the $ \gamma/Z $-penguin diagrams shown in Fig. 4(b). However, the resulting contributions to $ C_{9\ell}^ \rm{NP} $ from the γ- and Z-penguin diagrams are suppressed by the factors $ \alpha_e $ and $ 1-4s^2_W $, respectively. As a result, the dominant LFU-conserving contribution originates from the axial-vector part of the Z-boson couplings to fermions in the Z-penguin diagrams, with the final result given by

It can be seen from Eqs. (61) and (62) that for $ \kappa_{t, \nu} \sim \mathcal{O}(1) $ and $ m_S \sim \mathcal{O}(1)\, \mathrm{TeV} $, the LFU-violating coefficients $ C_{9\mu}^{\mathrm{NP}} = -C_{10\mu}^{\mathrm{NP}} $ numerically have the same order of magnitude as the LFU-conserving one $ C_{10\ell, Z}^ \rm{NP} $. Interestingly, this observation is also favored by the two-parameter fit for the $ R_{K^{(\ast)}} $ resolution [24],

obtained at the $ 1\sigma $ level. While a negative central value of the LFU-conserving coefficient $ C_{10\ell, Z}^{\rm{NP}} $ is, by itself, not helpful for explaining $ R_{K^{(\ast)}} $ anomalies, it can change the direction of the LFU-violating coefficients $ C_{9\mu}^{\mathrm{NP}} = -C_{10\mu}^{\mathrm{NP}} $ and, in particular, lift $ C_{9\mu}^{\mathrm{NP}}(=-C_{10\mu}^{\mathrm{NP}}) $ to a larger negative value compared to the one-parameter fits obtained in Refs. [18–28]. Translated to the parameter space in the $ t\nu $2HDM framework, this requires larger $ \kappa_{t, \nu} $ and/or a lighter scalar mass $ m_S $ to explain the $ R_{K^{(\ast)}} $ anomalies.

As can be inferred from the previous studies in Refs. [13–16], an $ \mathcal{O}(1) $ $ \kappa_\nu $ is generally required to explain $ R_{K^{(\ast)}} $ anomalies. Specific to the $ t\nu $2HDM framework, such a large coupling will readily help the right-handed neutrino $ \tilde \nu_R $ establish thermal equilibrium with the SM plasma via the Higgs doublet portal $ H_2 $. When the temperature T drops below the muon mass, the effective four-fermion interaction governing the right-handed neutrino annihilation rate, $ \Gamma_{2\tilde{\nu}\to 2\nu}\equiv \Gamma (\tilde{\nu}_{R} \bar{\tilde\nu}_{R} \to\nu_L \bar{\nu}_L) $, mediated by neutral scalars will determine the decoupling temperature $ T_{\tilde{\nu}, \rm dec} $ of $ \tilde{\nu}_R $. Because $ \tilde\nu_R $ is relativistic in the early Universe, its contribution to Hubble expansion can be parameterized by a shift in the effective neutrino number [31–35],

Here, $ N_{\tilde{\nu}}=1 $ denotes the number of thermalized right-handed neutrino species, and $ g_{\tilde{\nu}} =2 $ takes into account the antiparticle state of the right-handed Dirac neutrino, and $ g_{\tilde{\nu}}=1 $ for the right-handed Majorana neutrino. The effective d.o.f for the SM entropy density, $ g^s_\ast (T_{\tilde{\nu}, \rm dec}) $, is evaluated at the decoupling temperature $ T_{\tilde{\nu}, \rm dec} $ of $ \tilde{\nu}_R $, which can be estimated via the instantaneous decoupling condition $ \Gamma_{2\tilde{\nu}\to 2\nu}\simeq H(T_{\tilde{\nu}, \rm dec}) $, with the Hubble expansion rate given at the radiation-dominated epoch by

where the effective d.o.f for the energy density is taken approximately as $ g_*\equiv g_\ast^\rho \approx g_\ast^s $, and the Planck mass is given by $ M_{P} =1.22\times 10^{19}\, \mathrm{GeV} $.

The annihilation rate of the process $ \tilde{\nu}_{R} \bar{\tilde\nu}_{R} \to\nu_L \bar{\nu}_L $ has the structure

where $ |v_{\tilde{\nu}_{R}}-v_{\bar{\tilde\nu}_{R}}| $ is the relative velocity of the two incoming particles, $ g_{\tilde{\nu}}/2 $ is introduced to signify the symmetry factor due to the indistinguishability between the particle and antiparticle in the initial state, and $ n_{\tilde{\nu}} $ is the thermal particle number density of $ \tilde{\nu}_R $, given by

Here, note that the spin d.o.f of the right-handed neutrino $ \tilde{\nu}_R $ is equal to one for both the chiral Dirac and Majorana neutrinos. The thermal rate $ \langle\sigma_{2\tilde{\nu}\to 2\nu} |v_{\tilde{\nu}_{R}}-v_{\bar{\tilde\nu}_{R}}|\rangle $ in Eq. (66) is given by

where $ K_1(x) $ is the modified Bessel function of order one, and the phase-space factor is defined by

Within the $ t\nu $2HDM framework, the annihilation cross section $ \sigma_{2\tilde{\nu}\to 2\nu} $ is simply given by

where $ \sqrt{\hat s}=E_{\rm cm} $ is the center-of-mass energy. Finally, we obtain the annihilation rate of the process $ \tilde{\nu}_{R} \bar{\tilde\nu}_{R} \to\nu_L \bar{\nu}_L $,

which, together with the instantaneous decoupling condition $ \Gamma_{2\tilde{\nu}\to 2\nu}\simeq H(T_{\tilde{\nu}, \rm dec}) $ and Eq. (65), leads to the decoupling temperature

It should be mentioned that, to obtain the analytic thermal rate, as given by Eq. (68), we use the Boltzmann distribution $f^{\rm eq}_{\tilde{\nu}}={\rm e}^{-E/T}$. Because the dependence of the d.o.f $ g_*^s(T) $ on the decoupling temperature $ T_{\tilde{\nu}, \rm dec} $ is weak below the muon mass scale $ T<m_\mu $ [127], the effective neutrino number shift $ \Delta N_{\rm eff} $ will also have a weak dependence on $ T_{\tilde{\nu}, \rm dec} $. Therefore, the approximation of adopting the Boltzmann distribution is sufficient to estimate the scale of $ T_{\tilde{\nu}, \rm dec} $ from Eq. (72).

From Eq. (72), one can see that the decoupling temperature $ T_{\tilde{\nu}, \rm dec} $ will be solely determined by the free parameters $ \kappa_\nu $ and $ m_S $ after inserting the effective d.o.f $ g_*(T) $ as a function of the temperature [127]. This in turn implies that the effective d.o.f $ g_*^s(T_{\tilde{\nu}, \rm dec}) $ present in Eq. (64) and hence the effective neutrino number shift $ \Delta N_{\rm eff} $ are also determined by the two parameters $ \kappa_\nu $ and $ m_S $. However, given that the parameter $ \kappa_t $ is severely constrained by low-energy flavor physics (especially by the mass differences $ \Delta M_q $), we know that $ \kappa_{\nu} $ becomes the key parameter for the $ R_{K^{(\ast)}} $ resolution. Therefore, we can expect a potential correlation between the $ R_{K^{(\ast)}} $ resolution and the mitigation of $ H_0 $ tension achieved via the effective neutrino number shift given by Eq. (64) within the $ t\nu $2HDM framework proposed here.

Let us begin with the exploration of the NP parameter space allowed by $ R_{K^{(\ast)}} $ anomalies. By fixing the quasi-degenerate Higgs mass at $ m_S = 500 $, $ 700 $, $ 900 $, and $ 1200~{\mathrm{GeV}} $, we show in Fig. 5 the viable parameter regions in the $ (\kappa_\nu, \kappa_t) $ plane, under the perturbative unitarity bounds given in Eqs. (54) and (55). We also take into account all the relevant phenomenological constraints discussed in Sec. III. Explicitly, the regions above the various curves in Fig. 5 are already excluded by the branching ratio $ \mathcal{B}(\bar{B}\to X_s\gamma) $ (red), mass differences $ \Delta M_s $ (orange), and $ \Delta M_d $ (magenta), as well as by direct searches for the charged Higgs from ATLAS with $ 13\, \mathrm{TeV} $ (blue). In the upper two plots, we also show the correlation between $ \kappa_t $ and $ \kappa_\nu $ inferred from the CMS direct searches for neutral scalars, with the two benchmark points of the decay width over mass ratio, $ \Gamma_S/m_S = 10\% $ (black dashed) and $ \Gamma_S/m_S = 25\% $ (black solid). Because the LFU ratios of the di-lepton decays of $ Z/W $ gauge bosons do not impose any further significant constraints in the $ (\kappa_\nu, \kappa_t) $ plane under the perturbative unitarity bounds, they are not displayed in Fig. 5. Finally, the bands colored from dark to light green represent the regions allowed by the $ R_{K^{(\ast)}} $ resolution in the direction of $ C_{9\mu}^{\rm{NP}} = -C_{10\mu}^{\rm{NP}} $ at $ 1-3\sigma $, whereas the band in yellow denotes the $ 1\sigma $ region of $ C_{10\ell, Z}^{\rm NP} $, as given in Eq. (63).

Figure 5.  (color online) Viable parameter space in the $ (\kappa_\nu, \kappa_t) $ plane for the $ R_{K^{(\ast)}} $ resolution, with the Higgs mass fixed at $ m_S = 500 $, $ 700 $, $ 900 $, and $1200~ {\mathrm{GeV}}$. The green and yellow bands represent the regions allowed by the two-parameter fits [24] in the direction of $ C_{9\mu}^{\rm{NP}} = -C_{10\mu}^{\rm{NP}} $ and $ C_{10\ell, Z}^{\rm NP} $, as given by Eq. (63). We also take into account all the relevant constraints discussed in the previous two sections; see text for further details.

From Fig. 5, it can be readily seen that the most stringent bound in the quark sector originates from the mass differences $ \Delta M_q $ and, in particular, from $ \Delta M_d $, which in turn requires $ \kappa_\nu\simeq 3 $ for the $ R_{K^{(\ast)}} $ resolution. However, as mentioned in Sec. III.A.2, the $ \Delta M_q $ constraints may not be so conclusive owing to the large theoretical uncertainties. It should also be pointed out that if the $ \Delta M_q $ discrepancies observed in Sec. III.A.2 are confirmed in the future, we will have to introduce extra NP sources beyond the minimal $ t\nu $2HDM setup considered here. In such a special case, the constraints from $ \Delta M_s $ (orange) and $ \Delta M_d $ (magenta) may become irrelevant. However, for $ m_S = 500 $ and $m_S = 700~ {\mathrm{GeV}}$, the two black curves inferred from the CMS direct searches for neutral scalars should be interpreted as the maximal values of $ \kappa_t $ under the reference values of $ \kappa_\nu $. For instance, with $ \kappa_\nu\approx2 $, $ \kappa_t>0.97 $ will be excluded by the limits set by the CMS direct searches for the process $ pp\to S\to t \bar t $ [120]. Furthermore, for the benchmark point $ \Gamma_S/m_S = 10\% $ and $m_S=500~ \mathrm{GeV}$, the constraint on $ \kappa_t $ from the CMS neutral-scalar searches is tighter than that from the charged-Higgs bound set by the ATLAS collaboration with $ 13\, \mathrm{TeV} $ [123], whereas for $ \Gamma_S/m_S $ up to $ 25\% $ and $m_S=700~ \mathrm{GeV}$, the upper limit on $ \kappa_t $ is still determined by the CMS neutral-scalar searches. However, we must note that the CMS constraints are no longer applicable for $m_S > 750~ \mathrm{GeV}$ [120].

In the next subsection, we show that a large muon-related coupling $ \kappa_\nu $, as required by the $ R_{K^{(\ast)}} $ resolution, is necessary for generating a significant contribution to the $ \Delta N_{\rm{eff}} $ shift. In this respect, we conclude that the $ t\nu $2HDM framework provides us with an opportunity to correlate the $ R_{K^{(\ast)}} $ resolution with the mitigation of $ H_0 $ tension.

To visualize the potential correlation between the $ R_{K^{(\ast)}} $ resolution and the mitigation of $ H_0 $ tension achieved via the $ \Delta N_{\rm{eff}} $ shift, we start with Eq. (64), where the effective d.o.f $ g_*^s(T_{\tilde{\nu}, \rm dec}) $ is solely determined by the free parameters $ \kappa_\nu $ and $ m_S $ within our approximation. The LFU-violating Wilson coefficients $ C_{9\mu}^{\rm{NP}} = -C_{10\mu}^{\rm{NP}} $ in Eq. (61) can then be expressed in terms of $ \kappa_t $, $ m_S $, and $ \Delta N_{\rm{eff}} $. By fixing the scalar mass $ m_S $ and varying the parameter $ \kappa_t $ from zero up to the values allowed by the $ \Delta M_q $ constraints, we can finally obtain the numerical dependence of $ C_{9\mu}^{\rm{NP}} = -C_{10\mu}^{\rm{NP}} $ on $ \Delta N_{\rm{eff}} $, which is shown in Figs. 6(a) and 6(b) for $ m_S = 500 $ and $m_S = 1000~ {\mathrm{GeV}}$, respectively. The horizontal bands colored from dark to light green correspond to the global-fit results of $ C_{9\mu}^{\rm{NP}} = -C_{10\mu}^{\rm{NP}} $ at the $ 1- 3\sigma $ level, as given in Eq. (63). The blue (red) region corresponds to the viable parameter space in the Dirac (Majorana) neutrino case, in which the peak of $ C_{9\mu}^{\rm{NP}} = -C_{10\mu}^{\rm{NP}} $ corresponds to the upper limit on $ \kappa_{\nu} $, as required by the perturbative unitarity bound (see Eq. (54)).

Figure 6.  (color online) Induced ranges of $ C_{10\mu}^{\rm{NP}}=-C_{9\mu}^{\rm{NP}} $ for given values of $ \Delta N_{\rm{eff}} $, with the parameter $ \kappa_t $ varied from zero up to the values allowed by the $ \Delta M_q $ constraints, and the scalar mass fixed at $m_S = 500~{\mathrm{GeV}}$(a) and $m_S = 1000~ {\mathrm{GeV}}$(b). The regions in blue and red correspond to the Dirac and Majorana natures of the right-handed neutrino, respectively. The horizontal bands in green correspond to the global-fit results of $ C_{10\mu}^{\rm{NP}}(=-C_{9\mu}^{\rm{NP}}) $ at the $ 1- 3\sigma $ level, as given in Eq. (63).

From Fig. 6, it can be clearly seen that the resolution of $ R_{K^{(\ast)}} $ anomalies at the $ 1\sigma $ level requires the shift $ \Delta N_{\rm{eff}} \simeq 1.0 $ for a one-flavor right-handed Dirac neutrino and $ \Delta N_{\rm{eff}} \simeq 0.5 $ for a one-flavor right-handed Majorana neutrino, and any large or moderate deviation from the benchmark points of $ \Delta N_{\rm eff} $, although able to ease $ H_0 $ tension, cannot simultaneously resolve the $ R_{K^{(\ast)}} $ anomalies. In both cases, after fixing the Higgs mass, a large effective $ \Delta N_{\rm eff} $ is always required by a value $ \kappa_\nu \simeq 3 $ that almost coincides with the perturbative unitarity limit, and varying the value of $ \kappa_t $ only influences the Wilson coefficients $ C_{9\mu}^{\rm{NP}} = -C_{10\mu}^{\rm{NP}} $. Moreover, by comparing the two figures, we can find that increasing the Higgs mass will enlarge the viable space of $ \Delta N_{\rm eff} $ and simultaneously shrink the range of $ C_{9\mu}^{\rm{NP}} = -C_{10\mu}^{\rm{NP}} $. This indicates a preference of a lighter Higgs to address $ R_{K^{(\ast)}} $ anomalies while easing $ H_0 $ tension. In addition, such a difference between Dirac and Majorana neutrinos is expected owing to the different spinor natures of the neutrinos involved, that is, the Weyl spinor for the former and the Majorana spinor for the latter case. In terms of the favored $ \Delta N_{\rm{eff}} $ shifts inferred from Eqs. (5)–(8), we can then conclude that the eV-scale Majorana neutrino embedded in the $ t\nu $2HDM framework is able to address $ R_{K^{(\ast)}} $ anomalies and simultaneously ease $ H_0 $ tension, whereas the case with the one-flavor Dirac neutrino generates too large a $ \Delta N_{\rm{eff}} $ shift.

Finally, it can be demonstrated that, if more than one neutrino contributes to $ R_{K^{(\ast)}} $ anomalies via the box diagram shown in Fig. 4(a), the resulting $ \Delta N_{\rm{eff}} $ shift would be unacceptably large. As an example, let us consider the case in which there are two right-handed neutrinos with significant couplings to the muon lepton in Eq. (13). We must sum over the two flavors of $ \nu_R $ in Eq. (61), that is, $ |\kappa_\nu|^2\to |\kappa_{\nu, 1}|^2+|\kappa_{\nu, 2}|^2 $. Assuming that $ \kappa_{\nu, 1}\sim \kappa_{\nu, 2} $ and applying our previous finding of $ \kappa_\nu\simeq 3 $ for the $ R_{K^{(\ast)}} $ resolution, as inferred from Fig. 5, we can see that $ \kappa_{\nu, 1}\sim \kappa_{\nu, 2}\sim \mathcal{O}(3/\sqrt{2}) $ are required in this case. This means that the muon-related couplings $ \kappa_\nu $ can be reduced by a factor of $ 1/\sqrt{2} $ when explaining $ R_{K^{(\ast)}} $ anomalies with two right-handed neutrinos. However, such a parameter reduction cannot cause any significant increase in the decoupling temperature and, more importantly, any significant increase in the effective d.o.f $ g_*(T_{\tilde{\nu}, \rm dec}) $ in Eq. (72). This can be understood by the fact that enhancing $ T_{\tilde{\nu}, \rm dec} $ by a factor of $ 2^{2/3} $ can only increase $ g_*(T_{\tilde{\nu}, \rm dec}) $ by approximately $ 10\% $ [127]. From Eq. (64), we can see that the $ \Delta N_{\rm{eff}} $ shift is basically determined by the number of thermalized right-handed neutrino species. As a consequence, the correlation shown in Fig. 6 indicates that the $ R_{K^{(\ast)}} $ resolution with more than one thermalized right-handed neutrino would introduce a large $ \Delta N_{\rm{eff}} $ shift beyond that favored by Eqs. (5)–(8). This is the reason why we only introduce the one-flavor thermalized right-handed neutrino $ \tilde{\nu}_R $ into the early Universe within our framework.

The latest updated measurements from the LHCb [5] and SH0ES [55] collaborations have respectively strengthened the deviations of the LFU ratio $ R_{K} $ in rare semi-leptonic B-meson decays and the present-day $ H_0 $ parameter in the Universe. If confirmed with more precise experimental measurements and theoretical predictions, they could be tantalizing hints of NP beyond the SM. In this paper, we construct a simple flavor-specific 2HDM, dubbed the $ t\nu $2HDM, where significant NP effects arise only from the one-flavor right-handed neutrino and the top quark. Such a framework is only characterized by the three free parameters $ \kappa_t $, $ \kappa_\nu $, and $ m_S $ in the alignment limit with a quasi-degenerate Higgs mass spectrum.

The $ t\nu $2HDM can explain the long-standing $ R_{K^{(\ast)}} $ anomalies via one eV-scale right-handed Majorana neutrino or one right-handed Dirac neutrino under the most relevant constraints from low-energy flavor physics, the perturbative unitarity condition, and LHC direct searches. However, in contrast with the three-flavor right-handed neutrino scenarios considered in Refs. [13–16], an intriguing prediction resulting from the parameter space for the $ R_{K^{(*)}} $ resolution with a one-flavor scenario points toward a moderate shift in the effective neutrino number, $ \Delta N_{\rm eff}=N_{\rm eff}-N_{\rm eff}^{\rm SM} $, at the early BBN and late CMB epochs. It is then found that while the $ \Delta N_{\rm eff} $ shift predicted in the Dirac neutrino case is still at $ \Delta N_{\rm eff}\simeq 1.0 $ and hence disfavored by the CMB polarization measurements, the shift induced in the Majorana case is $ \Delta N_{\rm eff}\simeq 0.5 $, which coincides with the ranges from Eqs. (5)–(8) favored to ease the notorious $ H_0 $ tension [63–67]. There is also a potential correlation between $ R_{K^{(\ast)}} $ anomalies and the $ H_0 $ tension achieved via the $ \Delta N_{\rm eff} $ shift with the one-flavor eV-scale right-handed Majorana neutrino, and such a correlation can be tested in the future.

In conclusion, the $ t\nu $2HDM provides an interesting link between $ R_{K^{(*)}} $ anomalies and $ H_0 $ tension. In addition, a light right-handed Majorana neutrino embedded in the 2HDM infers a hierarchical Majorana neutrino pattern for the seesaw generation of the neutrino masses and, in particular, a nearly massless active neutrino.

As a final comment, we discuss direct searches of right-handed neutrinos. These right-handed neutrinos, which are also called heavy neutral leptons with masses above the eV scale, are often proposed to explain several puzzles of fundamental physics, a foremost example being neutrino oscillations. These hypothetical particles can be of Majorana or Dirac nature. The present generation of experiments usually focuses on the following three aspects: neutrino masses, oscillation parameters, and neutrinoless double beta decay [128, 129]. Future precise measurements of these parameters may come from many types of experiments, such as short-baseline, fixed-target, and collider experiments. With the upcoming precision era of neutrino physics, these terrestrial experiments are expected to determine the exact mixing pattern and flavor structures of heavy neutral leptons [128]. In addition, specific to the $ t\nu $2HDM, new interactions in the lepton sector can lead to the charged-Higgs decaying into right-handed neutrinos, $ H^+\to\mu^+\nu $. These right-handed neutrinos can therefore be searched for at the LHC in terms of SM-like Yukawa interactions with extended neutrinos. However, such processes have not yet been observed at the LHC, and only some phenomenological studies exist in literature [125]. We expect that right-handed neutrinos will be detected via the channel $ H^+\to\mu^+\nu $ in future experiments, and the free parameters related to heavy neutral leptons will be determined by forthcoming neutrino experiments.

We thank Biao-Feng Hou for providing us with the MadGraph5_aMC@NLO calculation and helpful discussions.