Combined explanations of B-physics anomalies: the sterile neutrino solution

In this paper we provide a combined explanation of charged- and neutral-current $B$-physics anomalies assuming the presence of a light sterile neutrino $N_R$ which contributes to the $B \to D^{(*)} \tau \nu$ processes. We focus in particular on two simplified models, where the mediator of the flavour anomalies is either a vector leptoquark $U_1^\mu \sim ({\bf 3}, {\bf 1}, 2/3)$ or a scalar leptoquark $S_1 \sim ({\bf \bar 3}, {\bf 1}, 1/3)$. We find that $U_1^\mu$ can successfully reproduce the required deviations from the Standard Model while being at the same time compatible with all other flavour and precision observables. The scalar leptoquark instead induces a tension between $B_s$ mixing and the neutral-current anomalies. For both states we present the limits and future projections from direct searches at the LHC finding that, while at present both models are perfectly allowed, all the parameter space will be tested with more luminosity. Finally, we study in detail the cosmological constraints on the sterile neutrino $N_R$ and the conditions under which it can be a candidate for dark matter.


Introduction
Various intriguing hints of New Physics (NP) have been reported in the last years in the form of lepton flavour universality (LFU) violations in semileptonic B decays. In particular the R(D ( * ) ) = B(B → D ( * ) τ ν)/B(B → D ( * ) ν) observable in b → cτ ν charged current transition, with = e, µ, has been measured by the BaBar [1,2], Belle [3][4][5] and LHCb [6][7][8] collaborations to be consistently above the Standard Model (SM) predictions. Once global fits are performed [9,10], the combined statistical significance of the anomaly is just above the ∼ 4σ level. Other deviations from the SM have been observed in the LFU ratios of neutral-current B decays, R(K ( * ) ) = B(B → K ( * ) µ + µ − )/B(B → K ( * ) e + e − ) [11,12]. Also in this case the overall significance is around 4σ. This discrepancy, if interpreted as due to some NP contribution in the b → s ¯ transition, is further corroborated by another deviation measured in the angular distributions of the process B → K * µ + µ − [13,14], for which however SM theoretical predictions are less under control.
Finding a combined explanation for both anomalies in terms of some Beyond the SM (BSM) physics faces various challenges. In particular, in the SM the b → cτ ν transition occurs at tree-level and an explanation of the R(D ( * ) ) anomaly generally requires NP close to TeV scale, for which several constraints from direct searches for new states at collider experiments as well as in precision electroweak measurements and other flavour observables can be stringent, see e.g. Refs. . On the other hand, the neutral current b → s + − transition occurs in the SM through a loop-induced process, thus hinting to a higher NP scale or smaller couplings responsible for the R(K ( * ) ) anomaly.
Concerning the R(D ( * ) ) observables, it has recently been proposed that the measured enhancement with respect to the SM prediction can also be obtained by adding a new right-handed fermion, singlet under the SM gauge group, hereafter dubbed N R [42,43] (see also [26,[44][45][46][47] for earlier related studies). Differently from other explanations where the NP contributions directly enhance the b → cτ ν τ transition, this solution allows to evade the stringent constraints arising from the SU (2) L doublet nature of the SM ν τ neutrino. In this case the B → D ( * ) τ ν decay rate becomes the sum of two non-interfering contributions: B(B → D ( * ) τ ν) = B(B → D ( * ) τ ν τ ) + B(B → D ( * ) τ N R ).
Several effective operators involving N R can be written at the B-meson mass scale. In order to ensure that the differential distributions in the B → D ( * ) τ N R process are compatible with the SM ones, as implicit in the global fits where the experimental acceptances are not assumed to be drastically modified by the presence of extra NP contributions, we assume that the sterile neutrino has a mass below ∼ O(100) MeV [43] and that the dominant contributions to the R(D ( * ) ) anomaly is given by a right-right vector operator Matching to the observed excess one finds [9] (Summer 2018 update [10]) where v ≈ 246 GeV is the vacuum expectation value of the SM Higgs field. This gives a NP scale required to fit the observed excess Λ/ √ c R D = (1.27 +0.09 −0.07 ) TeV .
Such a low NP scale strongly suggests that this operator could be generated by integrating out at tree-level some heavy mediator. There are only three possible new degrees of freedom which can do that: • a charged vector W µ ∼ (1, 1, +1), • a vector leptoquark U µ 1 ∼ (3, 1, +2/3), • a scalar leptoquark S 1 ∼ (3, 1, +1/3), where in parentheses we indicate their SU(3) c ×SU(2) L ×U(1) Y quantum numbers 1 . The case of the W µ has been recently studied in detail in Refs. [42,43]. In this work we focus on the two coulored leptoquark (LQ) models. Interestingly enough, both LQs can also contribute to the neutral-current b → sµ + µ − transition. In particular, the vector LQ U 1 contributes to that process at tree-level while the scalar S 1 only at one loop.
By considering the most general gauge invariant Lagrangians and assuming a specific flavour structure, we study in details the conditions under which the two LQ models can simultaneously explain both the R(D ( * ) ) and the R(K ( * ) ) measured values, taking into account all the relevant flavour and collider limits. Our findings show that the vector LQ provides a successful combined explanation of both anomalies, while being consistent with other low and high p T experiments. Instead, while the scalar LQ can address R(D ( * ) ), a combined explanation of also R(K ( * ) ) is in tension with bounds arising from B s −B s mixing. Also, by studying the present limits and future projections for collider searches, we find that the Large Hadron Collider (LHC) will be able to completely test both models already with ∼ 300 fb −1 of integrated luminosity.
For both models we then show that additional contributions to the mass of the active neutrinos generated by the operator responsible for reproducing the R(K ( * ) ) anomaly point to a specific extension of our framework, where neutrino masses are generated via the inverse see-saw mechanism [48][49][50]. We finally study the cosmological bounds on the right-handed neutrino N R and discuss the conditions under which it can be identified with a Dark Matter (DM) candidate. We show that an O(1) keV sterile neutrino can behave as DM only when the operators responsible for the explanation of the R(K ( * ) ) anomaly are turned off. In this case N R can reproduce the whole DM abundance observed in the Universe under the condition of additional entropy injection in the visible after the N R decoupling, while being compatible with bounds arising from the presence of extra degrees of freedom in the early Universe and from structure formations at small scales.
Very recently, while this work was already in the final stages of preparation, Ref.
[51] appeared on the arXiv which has some overlap with our paper. In particular [51] also studies explanations of R(D ( * ) ) anomalies with the two LQs considered here, as well as with other states which generate operators different than the right-right one, and studies the present LHC limits from LQs pair production. In this work we go beyond that analysis by studying in detail the possibility of a combined explanation with the b → s + − neutral-current anomalies, by studying also LHC constraints from off-shell exchange of LQs, which turn out to be very relevant, by discussing a possible scenario that can account for the generation of neutrino masses and by presenting a detailed study of the cosmological aspects of the sterile neutrino relevant for the anomalies.
The layout of the paper is as follows. In Sec. 2 we introduce the two LQ models with a right-handed neutrino and we describe their flavour structure and their implications for the relevant flavour observables. Limits arising from LHC searches are shown in Sec. 3, while possible model extensions that can account for the generation of neutrino masses are discussed in Sec. 4. Sec. 5 is dedicated to the discussion of the cosmological properties of N R . We finally conclude in Sec. 6.

Simplified models and flavour observables
In this Section we separately describe the interaction Lagrangians of the two candidate LQs, U 1 and S 1 in the presence of a right-handed SM singlet N R , assuming baryon and lepton number conservation. We work in the down-quark and charged-lepton mass basis, so that Integrating out the LQs at the tree-level one generates a set of dimension-six operators, charged-current anomalies can be addressed while at the same time being consistent with all other experimental constraints. Furthermore, we also consider the possibility of addressing with the same mediators the neutral-current R(K ( * ) ) anomalies.

Vector LQ U 1
The general interaction Lagrangian of the vector LQ U 1 ∼ (3, 1, +2/3) with SM fermions and a right-handed neutrino N R reads where g q,d are 3×3 matrices while g u is a 3-vector in flavour space. The integration of the U 1 state produces the seven dimension-six operators indicated in Tab. 1, where ξ = v 2 /(4m 2 U ). From these operators it is clear that this vector LQ can contribute to R(D ( * ) ) in several ways: ii) via the scalar operator O ledq proportionally to g d bτ g q b(s)τ ; iii) via the scalar operator O lN uq proportionally to g u cN g q bτ ; iv) via the vector RR operator O eN ud proportionally to g u cN g d bτ . The first three solutions involve a large coupling to third-generation left-handed quarks and leptons and have been studied widely in the literature [30,38,[52][53][54][55][56][57][58][59][60]. Such structures can potentially lead to some tension with Z boson couplings measurements, LFU tests in τ decays, and B s −B s mixing. To avoid these issues and since our goal is to study mediators contributing to R(D ( * ) ) mainly via the operator in Eq. (1), we set g q iτ ≈ 0 and focus instead on case iv). In order to explain both the R(D ( * ) ) and R(K ( * ) ) anomalies we assume the LQ couplings to fermions to have the following flavour structure: with g d bτ g u cN ∼ O(1), g q bµ , g q sµ 1. Note that one could potentially also add a coupling to the right-handed top, but since it does not contribute to the flavour anomalies we neglect it in the following.
By fitting the excess in the charged-current LFU ratios one obtains with this coupling structure With the couplings in Eq. (5), the vector LQ also contributes at the tree-level to b → sµ + µ − transitions via the two operators O 1,3 lq . By fitting the anomaly and matching to the standard weak Hamiltonian notation we get where we used the result of the global fit in [61] (see also [62][63][64][65][66][67][68]). This corresponds to The vector LQ, with the couplings required to fit the B-anomalies as detailed above, contributes also to other flavour and precision observables. While all constraints can be successfully satisfied, we list in the following the most relevant ones. The contribution to the B c → µN decay width and the corresponding limit [69] are given by where f Bc ≈ 0.43 GeV [70], m Bc ≈ 6.275 GeV and τ Bc ≈ 0.507 × 10 −12 s [71]. A chirallyenhanced contribution is also generated for the D s → µN decay, which is measured at a few percent level: where Λ cs The prediction for the lepton flavour violating (LFV) decay B s → τ µ from the (O ledq ) µτ bs operator is given by where f Bs ≈ 0.224 GeV [70], m Bs ≈ 5.37 GeV and τ Bs ≈ 1.51 × 10 −12 s [71]. The only weak constraint on this decay is the indirect one arising from the total lifetime measurements of the B s meson, but in the future this process could be directly looked for at Belle-II.
A contribution to B s −B s mixing is generated at the loop level and is proportional to (g q bµ (g q sµ ) * ) 2 , which makes it negligibly small given Eq. (9). These couplings also induce a tree-level contribution to b → cµν, which is constrained at the ∼ 1% level, however also the prediction for this observable is well below the experimental bound due to the small size of the couplings.
Finally we notice that at one loop the vector LQ generates also contributions to Z couplings to SM fermions, precisely measured at LEP-1. These effects can also be understood from the renormalisation group (RG) evolution of the operators in Tab. 1 from the scale m U down to the electroweak scale [72][73][74]. The relevant deviations in Z couplings are: 2 where the 95% confidence level (CL) limits have been taken from Ref. [75]. It is clear that the O(1) couplings required to address the R(D ( * ) ) anomalies do not induce any dangerous effects in these observables.
We conclude this section by stressing that the vector LQ U 1 with the coupling structure in Eq. (5) is able to successfully fit both charged-and neutral-current B-physics anomalies, while at the same time satisfying all other flavour and precision constraints with no tuning required. In Sec. 3 we show how this mediator can also pass all available limits from direct searches, but it should be observed with more data gathered at the LHC. Finally, in Sections 4 and 5 we show how the sterile neutrino N R can satisfy all constraints from both neutrino physics and cosmology.

Scalar LQ S 1
The general interaction Lagrangian for the scalar LQ S 1 ∼ (3, 1, +1/3) and a right-handed neutrino where λ q,u are 3 × 3 matrices while λ d is a 3-vector in flavour space and the supscript c denote the charge conjugation operator. The operators generated by integrating out this LQ are listed in Tab. 1. As for the vector LQ, also the scalar can contribute to R(D ( * ) ) in several ways, including via a large coupling to third generation left-handed quarks and leptons [20,24,27,31,34,35,38,[76][77][78][79][80][81][82], which however leads to tension with electroweak precision tests and B s −B s mixing [38,82]. We thus focus on the case where g q iτ 1 and where the leading contribution to b → cτ ν arises from the operator in Eq. (1).
Contrary to the vector LQ, the scalar one does not contribute to b → sµ + µ − at the tree-level. It does, however, at one loop [20] via box diagrams proportionally to the λ q sµ λ q bµ couplings. Our goal is thus to fit R(D ( * ) ) at tree-level via right-handed currents involving N R , while possibly fitting R(K ( * ) ) at one-loop with the corresponding couplings to left-handed fermions. In this spirit we require the following couplings to be non-vanishing: In the limit where one does not address R(K ( * ) ), i.e. λ q qµ ≈ 0, the only NP contribution to R(D ( * ) ) is given by the operator in Eq. (1): which further implies Thus with O(1) couplings also the scalar LQ should live at the TeV scale in order to explain the measured values of R(D ( * ) ) . In the more general case, the couplings in λ q in Eq. (15) induce also different contributions to R(D ( * ) ) which can be relevant since, as shown below, λ q bµ should be large if one aims to fit R(K ( * ) ): also induces a chirally enhanced contribution to the LFV process B c → τν µ L : The corresponding constraint makes the contribution of these couplings to R(D ( * ) ) in Eq. (18) subleading, simplifying then the contribution to charged-current anomalies to the expression in Eq. (16).
The couplings to quark and lepton doublets λ q qµ generate a b → cµν charged-current transition, which implies a violation of LFU in b → c ν processes which is however constrained at the percent level [83] δR µe b→c ≈ 0.03 Since, as shown below, in order to fit R(K ( * ) ) the coupling λ q * bµ has to be larger than 1, it is necessary to tune the parenthesis as This relation also suppresses the non-interfering contribution to the same observable from the (O 1,3 lN qd ) µcb operators. Note that this relation corresponds to aligning the coupling to t L µ L in the up-quark mass basis, so that the LQ has a much suppressed coupling to c L . The same couplings also induce a possibly large tree-level contribution to b → sν µ L ν µ L . The 95% CL limit on B(B → K * νν) fixes the upper bound where in the second step we used the condition in Eq. (22).
The loop contribution to B → K ( * ) µ + µ − is given by [20] ∆C µ Hence an O(1) λ q bµ coupling is needed to explain the R(K ( * ) ) anomaly. This is compatible with the constraint in Eq. (23) for m S 2 TeV.
As for the case of the vector LQ, the RG evolution of the effective operators down to the electroweak scale generates an effect in Z couplings. In this setup this is particularly relevant for the Zµµ one, due to the contribution proportional to y 2 t : At one loop, the couplings λ q bµ and λ q sµ also contribute to B s −B s mixing: Figure 1: 95% CL limits from flavour observables and Z couplings measurements on λ q bµ as a function of the scalar LQ mass. The green (yellow) region represents the parameter space which fits R(K ( * ) ) at 1σ (2σ).
It is clear that some tension is present with the value required to fit R(K ( * ) ), Eq. (25), for any value of m S . These limits are shown in Fig. 1. While the model is compatible with the experimental bounds on B s mixing within 2σ if the result from UTfit [84] is considered, the bound from Ref. [85] (see also Refs. [86,87]) excludes the R(K ( * ) ) solution, unless some other NP contribution to B s −B s mixing cancels the one from S 1 .

Collider searches
In Sec. 2 we have shown that in order to explain the observed value of R(D ( * ) ) both the vector and the scalar LQ should have a mass that, for O(1) value of the couplings, are around 1 TeV, thus implying the possibility of testing their existence in high-energy collider experiments. At the LHC LQs can be searched for in three main ways: i) they can be produced on-shell via QCD interactions; ii) they can be singly produced via their couplings to SM fermions; iii) they can be exchanged in the t-channel in qq scattering.
In this Section we will illustrate the main constraints arising from LHC searches on the two considered LQ models from both pair-production and off-shell exchange. Singleproduction modes, instead, while will be relevant in the future for large LQ masses, at present do not offer competitive bounds, see e.g. Ref. [88].

Pair-production
The interactions of Eq. (4) can be constrained in several ways by LHC searches. When produced on-shell and in pairs through QCD interactions, the LQs phenomenology is only dictated by the relative weight of their branching ratios. As we discussed in Sec. 2, the couplings g q sµ and g q bµ in Eq. (4) can give R(K ( * ) ) at tree-level, thus implying that they should be considerably smaller than g d bτ and g u cN , which are responsible for explaining R(D ( * ) ) also at tree-level, see Eq. (9) and Eq. (7). For this reason g q sµ and g q bµ can be neglected while studying the LHC phenomenology of the vector LQ. The relative rate of the dominant decay channels is thus set by the following ratio Regarding production, LQs can be copiously produced in pairs at the LHC through QCD interactions described by the following Lagrangian Here g s is the strong coupling constant, G a µν the gluon field strength tensor, T a the SU (3) c generators with a = 1, ..., 8 and κ is a dimensionless parameter that depends on the ultraviolet origin of the vector LQ. The choices κ = 0, 1 correspond to the minimal coupling case and the Yang-Mills case respectively. Barring the choice of κ, the cross-section only depends on the LQ mass 3 . For our analysis we compute the LQ pair production cross-section at LO in QCD with MadGraph5 aMC@NLO [89] through the implementation of the Lagrangian of Eq. (29) in Feynrules performed in [88] that has been made publicly available 4 .
The CMS collaboration has performed various analyses targeting pair produced LQs. In particular the analysis in [91], recently updated in [92], searched for a pair of LQs decaying into a 2b2τ final state setting a limit of ∼ 5 fb on the inclusive cross-section times the branching ratio for a LQ with a mass of 1 TeV. In the case of the 2c2N R final state, we can reinterpret the existing experiental limits on first and second generation squarks decaying into a light jet and a massless neutralino [93], for which the ATLAS collaboration provided the upper limits on the cross-sections for various squark masses on HEPData, which have then been used to compute the bounds as a function of the LQ mass 5 .
The bounds arising from LQs pair production searches are shown as green and blue shaded areas in Fig. 2 for κ = 0 (left panel) and 1 (right panel) in the m U − g d bτ plane. Here g u cN has been fixed to match the central value of R(D ( * ) ) according to Eq. (7). Also shown are the projections for a LHC integrated luminosity of 300 fb −1 , which have been obtained by rescaling the current limits on the cross section by the factor 300 fb −1 /L 0 , with L 0 the current luminosity of the considered analysis. All together we see that current direct searches are able to constrain vector LQs up to ∼ 1.3 TeV for κ = 0, and ∼ 1.8 TeV for κ = 1 when the dominant decay mode is into a 2c2N R final state, with slightly weaker limits in the case of an inclusive 2b2τ decay.

Off-shell exchange
From the Lagrangian of Eq. (4), and with the assumptions of Eq. (5), we see that other relevant constraints can arise fromcc → N R N R ,bb → τ τ andbc →τ N R processes which occur through the exchange of a t-channel LQ.
In particular,bc →τ N R directly tests the same interactions responsible for explaining the R(D ( * ) ) anomalies. The ATLAS collaboration published a search for high-mass resonances in the τ ν final state with 36 fb −1 of luminosity [95], which we can use to obtain limits in our model. To do this, we computed with MadGraph5 aMC@NLO the fiducial acceptance A and reconstruction efficiency in our model as a function of the threshold in the transverse mass m T , and used the model-independent bound on σ(pp → τ ν + X) × A × as a function of m T published in [95] to derive the constraints. We then rescale the expected limits on the cross section with the square root of the luminosity to derive the estimate for future projections. The present and future-projected limits in the m U vs. |g u cN g d bτ | plane derived in this way are shown in Fig. 3, together with the band showing the region which fits the R(D ( * ) ) anomaly. We notice that, while the present limits are still not sensitive enough to test the parameter space relevant for the anomalies, with 300 fb −1 most of the relevant space will be covered experimentally. Also, with more and more luminosity, this channel will put upper limits on the LQ mass (when imposing a successful fit of the R(D ( * ) ) anomaly). This complements the lower limits usually derived from pair-production searches.
The cc → N RNR channel gives rise to a fully invisible final state. In this case one can ask for the presence of an initial state radiation jet onto which one can trigger, thus obtaining a mono-jet signature. The CMS collaboration has performed this analysis for the case of a coloured scalar mediator connecting the SM visible sector with a dark matter candidate [96]. By assuming only couplings with the up type quarks, and fixing this coupling to one, they obtain a bound of 1350 GeV on the LQ mass. This corresponds to a parton level cross-section of ∼ 16 fb for p j T > 250 GeV, which we use as an upper limit on the monojet cross-section to set the limits on the vector LQ mass and couplings. For the bb → τ τ process, we impose the bound obtained in [97] and rescale it with the √ L factor in order to get the estimate for the projected sensitivity.
The current and projected constraints arising from the off-shell analyses are shown together with those from LQ pair production searches in Fig. 2. We observe that monojet and τ τ searches nicely complement direct searches for small and large g d bτ , respectively. Impressively, the off-shell search for τ N R , which exclude the region on the right of the contours, will completely close the parameter space already with 300 fb −1 of integrated luminosity, thus making this scenario falsifiable in the near future.

Pair-production
As for the vector case, also the interactions of the scalar LQ in Eq. (14) can be constrained in several ways. The on-shell production of a pair of scalar LQs is the dominant search channel at the LHC, which only depends on the LQ mass and branching ratios. 6 Since in Sec. 2 we showed that the couplings λ q sµ and λ q bµ of S 1 that are needed to fit R(K ( * ) ) might be incompatible (depending on the SM prediction considered) with the constraints arising from B s −B s mixing, we set them to zero for the forthcoming discussion. For LQ pair production searches the phenomenology of the scalar LQ is thus determined by the following ratio The CMS analysis [94] searches LQs decaying into the bbνν final state. This analysis can be directly applied to the case of the scalar LQ, given than the only difference with the decay mode targeted by the experimental analysis is the nature of the final state neutrino, which however does not strongly affect the kinematics of the event. For the 2c2τ final state no direct searches exist. The CMS analysis in [91], recently updated in [92], targets the bbτ + τ − decay mode and in principle cannot be applied to our scenario. We however observe that, for 100% branching ratios, the cross section in the analysis signal region (σ SR ) for the LQ → cτ or bτ cases is given by where c is the probability to mis-identify a c-jet as a b-jet, b is the b-jet tagging efficiency, [A × ] i is the acceptance for the considered final state and σ LQ Th. is the LQs pair production cross section. Since the kinematics of the event is not expected to change if a final state quark is a b-jet or a c-jet, the ratio of the number of events in the signal region for the case of the bτ and cτ final state is simply given by 7 i.e. the cross section is rescaled by a factor only dictated by the jet tagging efficiencies. In particular the upper limit on the cross section has to be divided by the factor in Eq. (32) which is smaller than 1. For concreteness we use the 70% b-tag efficiency working point of [91] from which we obtain c ∼ 20% [98]. The bounds arising from LQs pair production searches are shown as green and orange shaded areas in Fig. 4 (left) in the m S − λ d bN plane for the 2b2N R and 2c2τ final state respectively, where λ u cτ has been fixed to match the central value of R(D ( * ) ), see Eq. (17). We also again show the projections for a higher LHC integrated luminosity, namely 300 fb −1 . All together we see that current direct searches are able to constrain scalar LQs with a mass of ∼ 1 TeV when the dominant coupling is the one to bN while a weak constraints of ∼ 600 GeV can be set if the dominant coupling is the one to cτ , with these limits becoming ∼ 1.3 TeV and 1 TeV respectively for 300 fb −1 .

Off-shell exchange
Similarly to the vector LQ, also the scalar S 1 can be exchanged in t-channel in cb → τN R , bb → N RNR , and cc → τ τ processes. Also in this case the cb → τN R process directly tests the same couplings involved in the explanation of the R(D ( * ) ) anomalies. The experimental limits, and future projections, are obtained from the ATLAS analysis [95] in the same way as described for the vector LQ case. The derived limits in the m S −|λ d bN λ u cτ | plane, superimposed with the 68% and 95% CL intervals around the central values for R(D * ), are shown in the right panel of Fig. 4. Also in the scalar LQ case this search will put an upper limit on the LQ mass m S once the fit of the charged current flavour anomalies is imposed, and the high luminosity phase of the LHC with 3000 fb −1 of integrated luminosity will cover the whole relevant parameter space. The bb → N RNR final state can be constrained by monojet searches in an analogous way as done for the vector LQ. The excluded parameter space is shown as a purple region in the left panel of Fig. 4.
The limits on the cc → τ τ process can be obtained from the ones computed in [97] for bb → τ τ case (shown in the bottom panel of Fig. 6 of [97]) by taking into account the different parton luminosities for the two different initial state quarks. In particular, we approximate the R cb (ŝ) = L cc (ŝ)/L bb (ŝ) ≈ 2.5 ratio as constant and rescale the limit on the y bτ L coupling in [97] neglecting the interference of the signal with the SM background: limit(|λ u cτ |) ≈ limit(|y bτ L |)R 1/4 cb . The resulting excluded region is shown as a red region in the left panel of Fig. 4.
All together the current and projected constraints arising from these three analyses are shown together with the one arising from LQ pair production searches in the left panel of Fig. 4. We observe that τ τ searches nicely complement direct searches for small λ q bN while also in this case searches for τ N R , which again exclude the region on the right of the contours, will almost completely close the parameter space already with 300 fb −1 of integrated luminosity.

Neutrino masses and decays
The phenomenology of both the SM-like and sterile neutrino crucially depends on whether only the R(D ( * ) ) anomalies are addressed or if also the neutral-current ones are. This is particularly relevant for the vector LQ, since this state allows to explain both without any tension with flavour, precision, or collider constraints. For this reason in the following we Figure 5: Diagram responsible for generating a ν − N R Dirac mass term at one loop in the vector LQ model in case both charged-and neutral-current anomalies are addressed.
discuss both scenarios separately, stressing the main consequences for each of them.

Addressing only R(D ( * ) )
The operator responsible for reproducing the R(D ( * ) ) anomalies, Eq. (1), generates a Dirac mass term L ∼ m Dντ L N R + h.c. at two loops, where one can estimate [42,43] m D R(D ( * ) ) ∼ Such a small contribution to neutrino masses does not affect their phenomenology in a relevant way and therefore can be mostly neglected. In this scenario the leading decay mode for the heavy neutrino is N R → ν τ γ, which also arises at two loops from the same operator, with a rate (see Ref. [43] and references therein) which is much larger than the age of the Universe.

Addressing also R(K ( * ) )
If one wants to address also the neutral-current anomalies R(K ( * ) ) the situation becomes more complicated. In the following we focus on the model with the vector LQ, since it is the one which allows to do so without introducing tension with other observables. The chirality-flipping operators O lN uq induce a Dirac mass term between N R and ν µ at one loop, see Fig. 5, and with less suppression from light fermion masses: where we used the constraint in Eq. (10). Such large neutrino masses are of course incompatible with experiments. One possible solution is to finely tune these radiative contributions with the corresponding bare Dirac neutrino mass parameter, in order to get small masses. A more natural and elegant solution can instead be found by applying the inverse see-saw mechanism [48,49] (see also [50]). This was also employed recently in the context of the B-meson anomalies in Ref. [58]. In its simplest realisation, this mechanism consists in adding another sterile state 8S L with a small Majorana mass µ S and Dirac mass M R withÑ R . By defining n = (ν L ,Ñ c R ,S L ) t the mass Lagrangian L n = −1/2n M n n c can be written in terms of the following mass matrix Diagonalising the matrix, in the limit µ S m D < M R , the spectrum presents a light SM-like Majorana neutrino with mass and two heavy psuedo-Dirac neutrinos N R 1,2 with masses m N R ∼ (m D ) 2 + M 2 R and a splitting of order µ S . A small enough µ S can therefore control the smallness of the contribution to the light neutrinos without the need of any fine tuning. The mixing angle between the light neutrinos and the sterile one is given by where we used the (conservative) experimental bound of Ref. [99] for sterile neutrinos with masses m N R ∼ 10 MeV. Indeed, this limits puts a lower bound on the mass of the sterile neutrinos m N R 10 2 m D ∼ 1 MeV which is relevant for the cosmological analysis of the model.
In this case, the main decay modes of the sterile neutrino are N R → 3ν, ν µ e + e − via the mixing with ν µ and an off-shell Z boson exchange [43]: In this scenario N R decouples from the SM thermal bath at a temperature of ∼ 300 MeV (see next section), then becomes non relativistic and behaves like matter, comes to dominate the energy density after big bang nucleosynthesis (BBN), and decays into neutrinos and electrons before the epoch of matter radiation equality. This would generate a large contribution to the SM neutrino and electron energy densities before CMB, which is not cosmologically viable.
To avoid this problem N R should decay before BBN, which requires τ N R < 1s. Looking at the leading decay mode, Eq. (39), a simple way to achieve this is to increase both m N R ≈ M R and m D such that m N R 130 MeV and θ νµN ∼ 10 −3 (satisfying the limits from Ref. [99]). In this case a suitable short lifetime can be obtained. Such a mass of the sterile neutrino is close to the bound where it could potentially have an effect on the kinematics of B → D ( * ) τ ν. However a precise analysis of this scenario can only be performed with all details of the experimental analysis available. Interestingly, there are almost no constraints on θ νµN in the window of ∼ 30 − 40 MeV (roughly the mass difference between the charged pion and the muon, see for example [100]). This window provides an opportunity for a short enough lifetime of N R in this model. Future measurements by DUNE [101] and NA62 [100] will be able to test the scenarios with m N R 130 MeV and with m N R ∈ [30,40] MeV.
Another possibility is to add a mixing of N R with the τ neutrino, by adding a suitable Dirac mass term. In this case the lower limits on θ ντ N [99] are much weaker, allowing θ 2 ντ N 10 −3 for m N R ≈ 100 MeV and even larger ones for lighter masses. This allows to reduce even further the N R lifetime, while keeping the N R mass below the 100 MeV threshold.

Cosmology of N R
In this section we discuss cosmological bounds and opportunities in the presence of right handed neutrinos. As we saw in the previous section, if we only want to address the R(D ( * ) ) anomaly the right handed neutrino can be as light as 10 −3 eV and is cosmologically stable. Instead, if we also address the R(K ( * ) ) anomaly then it is much heavier and with a shorter lifetime. In particular we showed that it must decay before BBN in order to be a viable option. In this section we focus on the case where only R(D ( * ) ) is addressed and N R is cosmologically stable.

Relic density
Addressing only R(D ( * ) ), N R can be light and has a lifetime longer than the age of the universe. It therefore contributes to the DM relic density. Fitting the R(D ( * ) ) anomaly fixes the strength of the interaction of N R with the right handed b, c, τ . This in turn implies that N R was in thermal equilibrium in the early universe, and determines when it decoupled from the thermal bath. Solving the Boltzmann equation (see Appendix A) we find that N R freezes out at a temperature of ∼ 300 MeV, slightly above the QCD phase transition. Since m N R 100 MeV in order to explain R(D ( * ) ), it is relativistic at freeze-out. Its relic abundance today, assuming a lifetime longer than the age of the universe, is then [102,103] Here s 0 = 2891 cm −3 is the present entropy density and ρ c = 1.05 × 10 4 h 2 eV cm −3 the critical energy density [71]. We find a yield n s today which ranges between 8.3 × 10 −3 and 1.3 × 10 −2 , and correspondingly 9 g * S (T dec ) in the range between 35 and 60. For the sake of the estimates which follow, we take g * S (T dec ) = 50 as our reference value. We see that m N R ≈ 50 eV can account for the required amount of DM in the universe. However this is now a hot relic, and as such it is not consistent with structure formation. To make it comply with these bounds, we can simply lower its mass. For m N R eV, the right handed neutrino makes up less than 2% of the DM abundance and it is safely within the structure formation bound [104].

∆N eff
Such a light N R contributes to the number of effective relativistic species, N eff . The quantity ∆N eff is defined as the ratio of the energy density in dark radiation and that in one species of SM neutrino at the time of BBN, The ratio of the temperatures can be found using the total entropy conservation in the visible sector, just after the right-handed neutrino decoupled from the thermal bath [105]: Thus, from Eq. (41), we get which is within the experimental constraints [106]. We then conclude that a minimal model with a single right-handed neutrino N R lighter than an eV can explain the R(D ( * ) ) anomalies and evade all the relevant cosmological constraints. However N R can only be a small fraction of the DM in this case.

The dark matter option and entropy injection
We have shown that in the minimal scenario N R is a hot relic and can only constitute a small fraction of the observed DM energy density. It is interesting to explore the possibility of raising the N R mass to the keV range to make it a warm dark matter candidate. From Eq. (40) we see that m N R ∼ keV results in overclosure of the universe. We can then consider adding to the model a second heavier right-handed neutrino, χ R , whose decay produces enough entropy to dilute the abundance of N R [107,108] 10 . The dilution factor, defined as and the yield turns out to be slightly higher, see Appendix A. The value of g * S (T ) has a strong dependence on T when we are close to the QCD phase transition, as in this case. We use g * S (T dec ) = 50 in the estimates that follow. The reader should keep in mind that, while in the right ballpark, this number has some degree of uncertainty. 10 For a recent application of the entropy dilution in the models with right-handed neutrinos see [109][110][111].
modifies the relic density and ∆N eff as Note that we need D of order 20 if we want to push m N R to the keV range. In what follows we study if we can achieve such a dilution in a rather minimal setup. We assume that the heavier right-handed neutrino χ R , analogously to N R , is subject to the interaction We want χ R to decouple from the thermal bath at high temperature (but still below Λ χ , so the use of the effective interaction is justified), to come to dominate the energy density of the universe, then to decay and reheat the universe between 300 MeV (the decoupling temperature of N R ) and BBN. We discuss each step in turn. χ R decouples from the thermal bath when Γ = n σv H, with σ = λ 2 s 16πΛ 4 χ (here s is the centre of mass energy squared). Assuming χ R is relativistic at decoupling, we find and a yield Then, as the universe expands and the temperature decreases, χ R becomes non relativistic, and eventually dominates the energy density. It decays when Γ χ H χ , with the Hubble parameter and the decay rate into b, c, τ We find the reheat temperature, T after χ decay , assuming that the energy density of χ R is instantaneously converted into radiation at decay, This temperature must be above BBN, but below the N R decoupling temperature: 1 MeV < T after χ decay < 300 MeV .
The dilution factor can be expressed as [107,108] D = g * (T after χ decay )T 3 after χ decay g * (T before χ decay )T 3 before χ decay  Figure 6: Isocontours of the entropy injection D (black lines). The lower red area is excluded because the reheated temperature T after χ decay is below 1 MeV and the upper red area becase the reheated temperature is above right handed neutrino decoupling temperature. Red contours indicate T after χ decay in MeV. T χ , shown in the right axis, is the decoupling temperature of χ R from the interaction in Eq. (46). D is shown in Fig. 6 in the M χ vs. λ plane as black contours, where we see that the entropy injection factor can reach at most ∼ 100. It is instructive to trade the parameters M χ and λ for T χ and T after χ decay , using Eqs. (51), (47), (50). Then the expression for the D becomes which indicates that the maximal value can be achieved for the maximal decoupling temperature T χ and the minimal reheat temperature. As in our scenario we restrict to a decoupling temperature below the mediator mass, ∼ 1 TeV, the maximal entropy dilution that can be achieved is D max ∼ 100. If we consider a higher decoupling temperatur, the dilution factor does not improve. The reason is that above the mediator mass the cross section for the bc ↔ τ χ scattering process scales as 1/s, rather than s/Λ 4 χ . As a result the reaction rate n σv is linear in T , implying that the process is out of equilibrium at very high temperatures, and freezes in at lower temperatures. When the temperature drops below Λ χ we are back to the scenario we have studied above.
We are now in the position to assess whether such a dilution factor leads to a successful model. We see from Eq. (45) that we can raise m N R up to 5 keV and have N R contribute to the totality of dark matter energy density. However, for this mass the decay N R → νγ is too fast (see Eq. (34)) and excluded by X-ray measurements, which put a bound τ N R →νγ > 10 26−27 s in that region [112,113]. To avoid this bound, we should push m N R down to ∼ 1 keV. This is in some tension with constraints on warm dark matter from the Lyman-α forest (see for example [113]), which prefers a sterile neutrino heavier than 3-5 keV. However, due to the large entropy dilution, our N R is slightly colder and likely to comply with the Lyman-α bound also when m N R ∼ 1 keV. Further detailed studied are needed to confirm if this is the case.

Conclusions
The set of deviations from the Standard Model observed in various B-meson decays, from different experiments and in various different observables, is one of the most compelling experimental hints for BSM physics at the TeV scale ever obtained. Even more interesting is the possibility that all the observed deviations could be explained in a coherent manner by the same new physics. This has been the focus of a large effort from the theory community in recent years and several attempts have been put forward to achieve this goal. It became clear that this is not an easy task, in particular due to the fact that the large size of the required new physics effect to fit the R(D ( * ) ) anomalies generates tensions with either highp T searches or other flavour observables. In this spirit, it has become important to look for other possible solutions to the anomalies with different theoretical assumptions, which might help to evade the constraints. One such possibility is that the BSM operator contributing to R(D ( * ) ) does not involve a SM neutrino but a sterile right-handed neutrino N R . If the operator has a suitable right-right vector structure and the sterile neutrino is light enough, the the kinematics of the process remain SM-like and the solution is viable.
In this paper we study two possible tree-level mediators for such operator in a simplified model approach: the vector leptoquark U µ 1 and the scalar leptoquark S 1 . In the first part of the paper we explore the possibility that these mediators could generate both chargedand neutral-current B-physics anomalies. We find that the vector U µ 1 , which contributes to b → sµµ at the tree-level, provides a viable fit with no tension with any other flavour observable. The scalar S 1 , instead, contributes to the neutral-current process at one loop, thus requiring larger couplings to fit R(K ( * ) ). This generates a tension with the bound from B s -B s mixing which makes the combined solution of both class of anomalies from this mediator disfavoured. For both models we study the present constraints, and future projections, from direct searches at the LHC, including all relevant on-shell LQ pair-production modes as well as channels where the LQ is exchanged off-shell in the t-channel. We find that at present both scenarios are viable, but already with 300 fb −1 of luminosity LHC will test almost all the viable parameter space. In particular, the search in the τ ν final state, which directly test the interactions relevant for the R(D ( * ) ) anomalies, puts upper limits on the LQ mass and in the future will completely cover the region which fits the anomalies.
In the second part of the paper we study the phenomenology of the sterile neutrino N R . This depends crucially on whether or not both classes of anomalies are addressed or only the charged-current ones are. In the former case a Dirac mass term with the muon neutrino is generated at one loop with a size of tens of keV. In order to keep the SM neutrinos light it is possible to employ the inverse see-saw mechanism, by introducing another sterile neutrino with a small Majorana mass and a large Dirac mass with N R . The outcome of this is that the SM-neutrinos are light but the sterile ones are above 10 MeV. The mixing between the muon and sterile neutrino induces a fast decay of N R , rendering it unstable cosmologically. To avoid issues with the thermal history of the Universe it should decay before BBN, which requires its mass to be ∼ 100 MeV.
If instead only the R(D ( * ) ) anomalies are addressed the picture changes completely. In this case a Dirac mass term with the tau neutrino is generated at two loops and it is small enough to not have any impact in neutrino phenomenology. The main decay of N R in this case is into ν τ γ and arises at two loops as well, with a lifetime much longer than the age of the Universe. In order not to overclose the Universe energy density its mass should be below ∼ 50 eV, which makes it a hot relic. The constraints on the allowed amount of hot dark matter impose an upper limit on its contribution to the present dark matter density, which translates into an upper bound for the mass m N R eV. If the sterile neutrino is to constitute the whole dark matter, an entropy injection at late times is necessary in order to dilute its abundance. This can be obtained, for example, by adding another heavy sterile neutrino which decays into SM particles after N R decouples. In this case we find that N R could be a warm dark matter candidate with a mass ∼ (few keV). This option is highly constrained by current cosmological and astrophysical observations. While our model seems to have a small region of viable parameter space, a conclusive statement requires further detailed studies.
To conclude, the U µ 1 model presented in this paper allows to fit both charged-and neutralcurrent anomalies with no tension at all with present low-and high-p T bounds. The sterile neutrino in this case is cosmologically unstable, decaying before BBN happens. In case one aims at only solving the R(D ( * ) ) anomalies, instead, the neutrino is stable and if it is light enough it satisfies all cosmological constraints. With some additions to the model, in particular a mechanism for entropy injection after it decouples, it can also be a candidate for dark matter at the keV scale.
involved. We use the following conventions, inspired by Ref. [114], z ≡ m b T , Y i = n i s e , n eq i,rel = g i T 3 π 2 , s e = g * S 2π 2 45 where we are using the Maxwell-Boltzmann statistics for simplicity. Here g i is the number of internal degrees of freedom of the particle (2 for a Weyl fermion), K 1 is a Bessel function, and with s the centre of mass energy squared. Depending on the mediator in the UV completion, one will also have effective operators which introduce either the N R N R ↔ cc, N R N R ↔ bb, N R N R ↔ τ τ scattering processes. Particularly important is the N R N R ↔ cc since charm is lighter than τ and b quark and is less Boltzmann suppressed, keeping N R in thermal equilibrium for a little longer. As a result, the effect of including the N R N R ↔ cc process is to slightly delay the freeze-out of N R . To take it into account we can add the term to the right hand side of Eq. (55). We show in Fig. 7 how Y N evolves as a function of z. We fix the interaction strength to Λ/ √ c R D = 1.27 TeV, which is the value which fits the R(D ( * ) ) anomaly, see Eq. (3). When the only processes are those in Eq. (55), we find the freeze-out temperature T FO,N 350 MeV , the final yield Y N,0 = 8.3 × 10 −3 , and When we include also the processes of Eq. (62) we find Note that these values of g * S should be taken with a grain of salt, as we are close to the QCD phase transition and g * S has a strong dependence on the temperature in this range.
The quoted values are meant as a ballpark which we use for the estimates in this paper.  (55), while the dashed line includes also the contribution from Eq. (62). We see that freeze-out occurs at z FO,N 12 (solid line) for processes involving only one N R , at z FO,N N 17 (dashed line) when we also include the process N R N R ↔ cc [see Eq. (62)].

Analytic estimates
We can check analytically the numerical result obtained above. Let's consider only the process N R τ ↔ bc. The equation of Boltzmann above is easily manipulated into the familiar formṅ N + 3Hn N = (−n N n eq τ + n eq N n eq τ ) σv , with σv ≡ γ(N τ → bc) n eq N n eq τ .
Written in terms of s (centre of mass energy squared) the rate density is with We know that at T = √ s min ∼ 5 GeV, for interactions not so much weaker than the weak force (that is for Λ in the TeV ballpark), N R is in thermal equilibrium. Thus, to make analytic progress, we can take the limit T √ s. In this limit