The role of right-handed neutrinos in $b \to c \tau \bar{\nu}$ anomalies

Motivated by the persistent anomalies reported in the $b\to c\tau\bar{\nu}$ data, we perform a general model-independent analysis of these transitions, in the presence of light right-handed neutrinos. We adopt an effective field theory approach and write a low-energy effective Hamiltonian, including all possible dimension-six operators. The corresponding Wilson coefficients are determined through a numerical fit to all available experimental data. In order to work with a manageable set of free parameters, we define eleven well-motivated scenarios, characterized by the different types of new physics that could mediate these transitions, and analyse which options seem to be preferred by the current measurements. The data exhibit a clear preference for new-physics contributions, and good fits to the data are obtained in several cases. However, the current measurement of the longitudinal $D^*$ polarization in $B\to D^*\tau \bar\nu$ cannot be easily accommodated within its experimental $1\sigma$ range. A general analysis of the three-body $B\to D \tau \bar\nu$ and four-body $B\to D^*(\to D\pi)\tau \bar\nu$ angular distributions is also presented. The accessible angular observables are studied in order to assess their sensitivity to the different new physics scenarios. Experimental information on these distributions would help to disentangle the dynamical origin of the current anomalies.

These experimental facts suggest a surprisingly large violation of lepton-flavour universality, and have triggered a large number of detailed phenomenological studies trying to determine the most plausible NP explanation. A quite complete list of relevant references can be found in Ref. [31], where an exhaustive analysis of all available data has been accomplished with a model-independent effective field theory (EFT) approach, assuming only the SM particle content and symmetries in order to define the basis of allowed low-energy operators. A global fit to all data, with a good statistical quality, has been obtained in terms of the four possible Wilson coefficients; however, the fit does not allow to clearly identify a potential mediator of the underlying NP interaction [31]. Moreover, the experimental value ofF D * L cannot be accommodated within 1σ [31].
Light right-handed neutrinos (RHNs) have been suggested [32][33][34][35][36][37][38][39][40][41][42] as a possibility to evade the current phenomenological constraints on the EFT operators containing left-handed neutrino (LHN) fields. Sterile neutrinos are singlets under the SM gauge group and, therefore, their properties are not linked to any charged electroweak partners. Moreover, the existing limits from the neutrino sector do not constrain significantly the scale of ν R operators beyond what is probed in b → cτν transitions. In order not to disrupt the measured B → D ( * ) τν invariant-mass distributions [4,6], one just needs to assume the ν R fields to be light, m ν R O(100) MeV, which also helps to avoid other cosmological and astrophysical limits. Neglecting neutrino masses, there is no interference between the two neutrino chiralities, and the decay probability becomes an incoherent sum of ν L and ν R contributions: B(b → cτν) = B(b → cτν L ) + B(b → cτν R ). Therefore, it is not difficult to increase the predicted rates towards the experimentally favoured range. However, a large ν R contribution requires the corresponding Wilson coefficients to be large, of the order of the SM ν L interaction, because the rates are quadratic in the ν R transition amplitude.
Previous works considering RHNs in B → D ( * ) τν decays [32][33][34][35][36][37][38][39][40][41][42] have focused on reproducing the integrated rates, most of them within particular scenarios of NP. All phenomenological analyses need to rely on the underlying assumption that the differential decay distributions, and hence the experimental acceptances, are not significantly modified by the NP contributions. While this assumption is unavoidable, in the absence of direct access to the data, none of the previous studies have included the measured q 2 distributions in their fits. This shape information has been shown to play an important role, discarding many proposed solutions with ν L fields [31,[43][44][45][46], and could be expected to be even more relevant for those solutions based on RHNs, since they induce distortions in the rates that are quadratic in NP contributions.
We aim to improve the situation in this paper, by extending the EFT analysis of Ref. [31] to a basis of dimension-six operators that includes light RHNs. In our fit procedure, we consider all observables measured for B → D ( * ) τν decays until date; including the data for binned differential distributions with respect to the lepton-neutrino invariant-mass squared, the D * longitudinal polarization fractionF D * L , the lepton polarization asymmetryP D * τ and the experimental results for R D ( * ) . The last ratios have been recently altered, reducing the tension with the SM and making a fresh re-analysis necessary. We also study the differential three-body B → Dτν decay distribution and derive the four-body angular distribution of the B → D * (→ Dπ)τν decay for the most general dimension-six Hamiltonian. By identifying the possible high-scale NP mediators which can generate the operators involving RHNs, we predict several angular observables that can be tested at the experiment.
The rest of the paper is organized as follows. In Section 2 the most general effective Hamiltonian for our analysis is described, and expressions of the relevant observables are written in terms of the Wilson coefficients. In Section 3 the experimental status of the b → c transitions is interpreted from an EFT approach, by looking at the effect that individual Wilson coefficients may produce in the relevant observables. In addition, all possible NP mediators that can effectively generate a b → cτν R transition, and the corresponding Wilson coefficients that will arise at low energies after their integration, are listed. In Section 4 the results of our fits are presented and discussed. We consider different scenarios, originated by the integration of the relevant NP mediators, and compare their fitted results with the SM case. Section 5 contains the predicted angular coefficients of the B → Dτν and B → D * (→ Dπ)τν distributions for the best fit scenarios, including the forward-backward asymmetries A D ( * ) F B , the τ polarization asymmetries P D ( * ) τ , and the integrated longitudinal polarization fraction F D * L . Finally, conclusions are exposed in Section 6. Many technical details, such as hadronic matrix elements, FFs, and the full set of relevant helicity amplitudes, are compiled in several appendices.
2 Theoretical framework and observables 2.1 Effective field theory with the ten four-fermion operators: which are invariant under SU (3) C ⊗ U (1) em . Tensor operators with different lepton and quark chiralities vanish identically. 1 The SM charged-current contribution to O V LL , from a W µ exchange, has been explicitly added to Eq. (5), so that C X AB = 0 in the SM. Any non-zero contribution to these Wilson coefficients is then a manifestation of NP beyond the SM. We are assuming that NP contributions are only present in operators involving charged leptons of the third generation. This is well justified, since potential NP effects have been shown to be negligible in b → c ν transitions [18].
In the subsequent sections, we present the analytic expressions of observables constructed from different decay modes containing the b → cτν quark-level transition. The presence of RHN only modifies the leptonic currents; therefore, the decay amplitudes can be given in terms of the same hadronic FFs used for ν L operators, that were already listed in Ref. [31]. For completeness, we compile them again in the Appendix A. We also list in Appendix B the whole set of relevant helicity amplitudes, which we have obtained following the standard helicity formalism for semileptonic B decays [47][48][49]. We have checked that our expressions reproduce the available results for LHNs, and that they satisfy the correct parity relations between the ν L and ν R transition amplitudes.

Observables
The relevant observables for our analysis can be classified into those involving B → D and B → D * semileptonic transitions. We also consider the leptonic decay B c → τν since it can constrain certain Wilson coefficients.

B → Dτν
The differential distribution of the decay B → Dτν can be written as where q 2 ≡ (p τ + pν) 2 , θ τ is the polar angle of the τ momentum in the rest frame of the τν pair, with respect to the z-axis defined by the momentum of the D meson in the B rest frame, and we have introduced the shorthand notation for the Källen function The coefficient functions of the different angular dependences are given by 1 This is a direct consequence of the Dirac-algebra identity γ5σ µν = − i 2 ε µναβ σ αβ , which implies σµν ⊗σ µν γ5 = σµν γ5⊗σ µν and σµν γ5 ⊗ σ µν γ5 = σµν ⊗ σ µν . We use the convention ε 0123 = −ε0123 = −1.
The hadronic helicity amplitudes H s V,0 , H s V,t , H s S and H s T are functions of q 2 , and their explicit expressions are given in the Appendix A. The LHN contributions to Eq. (9) are in full agreement with Ref. [50]. Notice that the vector and scalar Wilson coefficients only appear in the combinations C V LX +C V RX and C S LX +C S RX , regardless of the neutrino chirality X.
Integrating over cos θ τ , one obtains [38] dΓ The linear term in the θ τ distribution given in Eq. (7) can be accessed via the forward-backward asymmetry, traditionally defined as and the τ polarization asymmetry can be constructed as where the decomposition of the amplitude in τ helicity states is given in the Appendix B.3.

B → D * τν
The vector meson D * in the final state provides additional observables compared to the previous case. The angular analysis of a four-body final state, namely B → D * (→ D π)τν, further allows us to construct a multitude of observables that can be extracted from data. The differential decay distribution of the transition process on the mass shell, can be expressed in the form [50]: I s 1 sin 2 θ D + I c 1 cos 2 θ D + I s 2 sin 2 θ D + I c 2 cos 2 θ D cos 2θ τ + (I 3 cos 2φ + I 9 sin 2φ) sin 2 θ D sin 2 θ τ + (I 4 cos φ + I 8 sin φ) sin 2θ D sin 2θ τ + (I 5 cos φ + I 7 sin φ) sin 2θ D sin θ τ + I s 6 sin 2 θ D + I c 6 cos 2 θ D cos θ τ .
In addition to the lepton-pair invariant-mass squared q 2 = (p τ + pν) 2 , we use as kinematic variables the three angles φ, θ τ and θ D , which are defined as follows. Taking as positive zaxis the direction of the D * momentum in the B rest frame, θ τ and θ D are the polar angles of the τ and the final D meson in the τ ν and Dπ rest frames, respectively. The azimuth φ is the angle between the decay planes formed by τ ν and Dπ. See Fig. 1 for a visual representation of these kinematical variables. Measuring this four-dimensional distribution is obviously a major experimental challenge, since the subsequent τ decay involves one (τ → ν τ + hadrons) or two (τ → ν τ ν ) additional neutrinos, making difficult to reconstruct the τ direction. Some information can be recovered by measuring the distribution of the secondary τ decay [32], but we refrain to enter here into this type of technical (but important) details.
The angular coefficients I i 's are functions of q 2 that encode both short-and long-distance physics contributions. They can be written in terms of the hadronic FFs given in Appendix A. Using the global normalization where B(D * → D π) is the branching fraction of the D * decay into D π states described in Appendix C.
The expressions for the angular coefficients are: In the above expressions, the A L,R λ denote the transversity amplitudes, which are the projections of the total decay amplitude into the explicit polarization basis. The contribution of the RHN transitions to the angular coefficients is equivalent to the LHN ones, i.e. (L → R), up to a sign that depends on the relation between right-handed and left-handed leptonic transversity amplitudes. In the SM, the decay B → D * τν can be described by a total of four transversity amplitudes that correspond to one longitudinal (A 0 ) and two transverse (A ⊥, ) directions, and a time-like component (A t ) for the virtual vector boson decaying into the τν pair. However, with the inclusion of RHNs, we must distinguish the left and right chiralities of the leptonic current; thus, we get in total eight amplitudes: A L,R 0,⊥, ,t . Now, in presence of the NP operators given in Eq. (6), the (axial)vector contributions can be incorporated in the above mentioned eight transversity amplitudes, modified by the presence of the new Wilson coefficients. Nevertheless, the (pseudo)scalar and tensor operators induce eight further amplitudes (four for each neutrino chirality): two (pseudo)scalar amplitudes A L,R P and six tensor transversities A L,R T 0,T ⊥,T . Thus, with the most general dimension-six Hamiltonian in Eq. (5), the decay B → D * (→ D π)τν can be described by a total of sixteen tranversity amplitudes. Their explicit dependence on the hadronic helicity amplitudes, compiled in Appendix A, and the Wilson coefficients is listed below, where the t and the P amplitudes arise in the B → D * observables combined as With these definitions, the left-handed contributions to the angular coefficients in Eq. (16) are in agreement with Ref. [50]. Performing the angular integrations in Eq. (14), one easily obtains the differential distribution with respect to q 2 , given by which written explicitly in terms of the different Wilson coefficients takes the following form: Differential distributions with respect to a single angle, which can be obtained by integrating two angles at a time, are also of special interest. These are In the following we define several observables constructed from the coefficients of various angular dependences. The distribution with respect to cos θ D in Eq. (22) provides the longitudinal and transverse polarization fractions for the D * meson, defined as [50] which satisfy that F D * L + F D * T = 1. Notice that these quantities are functions of q 2 . The Belle measurement mentioned before in Section 1, named asF D * L , refers to the q 2 -integrated polarization. We define the q 2integrated observables as follows,Ō where Γ is the total decay width and the q 2 dependence of the observables has been written explicitly. The angular coefficients I 3 and I 9 can simply be extracted by measuring the terms proportional to cos 2φ and sin 2φ in Eq. (23), respectively. Furthermore, we define several asymmetries starting with the well-known forward-backward asymmetry, defined as The coefficients I 4 and I 5 in Eq. (14) can be extracted with the two angular asymmetries: One can further define the following two observables, which are non-vanishing only if NP induces a complex contribution to the amplitude. This holds true for the coefficient A 9 as well. These asymmetries are simply related to the angular coefficients in (14): Finally, the total branching ratio can be decomposed in terms of the τ polarization, giving rise to another observable: the lepton polarization asymmetry, defined as

B c → τν
Another interesting but yet not observed decay is B c → τν for which the branching ratio can be written as The first line contains the usual contribution from LHNs [31], while the ν R contribution in the second line is readily obtained with a parity transformation.

Interpreting the anomalies
This section is devoted to study the origin of the observed experimental deviations from the SM predictions. We show from a theoretical perspective the implications of new physics in the observables involving b → c transitions and discuss the possible ultraviolet (UV) scenarios that could give rise to such anomalies in the context of b → c processes involving both left-and right-handed neutrinos.

Fit-independent results
The Wilson coefficients introduced in Eq. (5) encode all NP contributions that can enter in b → c transitions at dimension-six operator level, also in the presence of sterile light RHNs. Therefore, the landscape of possibilities generating the anomalies can be classified by the impact of these ten parameters on the measurable observables. To get a general idea about the sensitivity to the different Wilson coefficients, we quote the numerical expressions of several observables that have already been measured. These expressions have been obtained setting the FFs at their central values and, therefore, ignoring the uncertainties and correlations among the different numerical factors. The complete analytical expressions, with a proper account of hadronic uncertainties, will be used instead in the data fits that we will present in Section 4. The observables R D and R D * are normalized to their SM predictions: For the q 2 -integrated polarization observablesP D * τ andF D * L , we show their numerical values multiplied by R D * :  Different Wilson coefficients could help to reproduce the measured values of R D and R D * . However, the scalar coefficients would need to take values that are already excluded by B(B c → τν), leaving vector and axial-vector contributions as the preferred options to fit the experimental results. The large uncertainties in theP D * τ measurement make almost any shift in the Wilson coefficients to be in agreement with the experimental value, being the only exceptions large shifts in the vector Wilson coefficients C V LR,RR and a positive increment of C T LL . Looking at the dependence of these observables on the RHN contributions, one observes that all of them are symmetric under the exchanges C V LR ↔ C V RR and C S LR ↔ C S RR . In particular,F D * L is insensitive to any combination of ν R vector contributions because the dependence on their corresponding Wilson coefficients exactly cancels, since it is defined as a ratio as Eq. (24) shows. This does not hold true for C V RL , since there is an interference between this NP operator and the SM contribution.
It is particularly challenging to reproduce the experimental value ofF D * L , regardless of the type of NP contribution; the ±1σ band cannot be reached varying any of the Wilson coefficients individually. Negative non-zero values of C V RL can only slightly increase the predicted longitudinal D * polarization, while the changes induced by the tensor Wilson coefficients go in the opposite direction of the experimental value, decreasing the SM predictions. The only contributions that would help are the scalar ones, but for values of their Wilson coefficients that are already excluded by the constraint B(B c → τν) < 30%.

UV Physics
Once the impact of individual Wilson coefficients in B → D ( * ) τν observables is understood, the following step is to extend the analysis to the combined effect of several coefficients that are present in these transitions simultaneously. The most general EFT Hamiltonian in Eq. (5) includes 10 Wilson coefficients, which in general can be complex. Even assuming them to be real, a 10-parameter fit would become unstable. Moreover, its interpretation in terms of NP mediators and UV completions might be unrealistic. Instead, we consider particular cases, described in Section 4.2. Most of them are motivated from the "simplified model" scenarios. In this context, "simplified" refers to a single new mediator particle that can be integrated out to contribute to one or more of the effective operators entering into the b → cτν transitions. As the main purpose of this work is to explore the effect of light RHNs, we single out those mediators that can contribute to the b → c transitions and involve a gauge-singlet RHN.
These NP fields can be classified into scalars, vector bosons and leptoquarks, as listed in Table 1. Since in most cases both right-and left-handed neutrino operators are generated simultaneously after a given mediator is integrated out, we will explore both the effect of considering only the right-handed contributions as well as the scenarios in which the full set of operators is generated. Unlike in previous references discussing the role of RHNs in b → c anomalies [35,38], we also include a fit to the Wilson coefficients that will appear if NP is mediated through the leptoquarkṼ 2 ∼ (3, 2, −1/6).

Fit Results
Under the assumption that NP enters only in the third generation of leptons and that Wilson coefficients are real, we have performed fits in different scenarios of the most general dimension-six Hamiltonian, taking into account all experimental data available nowadays. We start by listing the inputs used in the fit, and then we describe the motivated scenarios, based on the previous section, that we are considering. Finally, the results obtained by performing global fits in each of the scenarios are interpreted. Table 1:

Numerical input of the fit
For our fits we will use the most recent world-average values of R D and R D * from Ref. [13], including a correlation of −0.38 between them. The value of the q 2 -integrated τ polarization,P D * τ , measured recently by Belle [7] and the longitudinal D * polarization,F D * L , measured by BaBar [30] are also taken into account. Finally we consider the q 2 distributions of the D and D * meson [4,6], summarized in Table 9 of Ref. [31]. The different experimental inputs used in the fits are collected in Table 2. The upper bound for the leptonic decay rate B(B c → τν) is taken to be either 30% or 10%. The first limit is derived from the B c lifetime [46,51,52], while a stronger bound of 10% is obtained from the LEP data at the Z peak [53]. 2 In our analyses the stronger 10% limit is first assumed in the fit and, in those cases where the 10% bound is saturated the fit is repeated by relaxing it to 30%. As Eq. (32) shows, the B c → τν limit constrains the splittings between the C V LL(RR) and C V RL(LR) and, specially, between the C S RL(LR) and C S LL(RR) Wilson coefficients.
For the FFs, we follow the same approach as in Ref. [31]. Using heavy quark effective field theory [55,56], we adopt the Boyd, Grinstein and Lebed (BGL) parametrization [57][58][59], including corrections of order α s , Λ QCD /m b,c [15] and partly Λ 2 QCD /m 2 c [18]. We also include the cubic term in the expansion of the leading Isgur-Wise function in powers of the conformal-mapping variable z [60,61]. The inputs of the FFs, listed in Table 1 of Ref. [31], haven been obtained from lattice quantum chromodynamics (QCD) [62][63][64][65], light-cone sum rules [66] and QCD sum rules [67][68][69], without making use of experimental data. See Ref. [18] for more details on the FF parameters. As in Ref. [31], we allow the FF parameters to fluctuate around these input values, which are considered as pseudo-observables with their corresponding χ 2 FF taken into account in the fits. Since this theoretical χ 2 gives a very small contribution to the total χ 2 of the fits, we will not discuss it again and refer to Ref. [31] for additional technical details.

Scenarios and fit results
As previously mentioned, by adding RHN, the set of operators increases from 5 to 10. The large number of free parameters makes difficult to perform a global fit to the full basis of operators. Instead, we will work in different motivated scenarios that arise by integrating out a single NP mediator and, therefore, contribute to small subsets of operators at the m b scale. Possible candidates, their quantum numbers and the operators generated once the given mediator is integrated out are listed in Table 1. The last two columns show the operators involving left-handed and right-handed neutrinos. Following previous works, we consider scenarios that only take into account the contributions from RHN operators, labelling them with the letter "a" [35], while "b" scenarios also contain the LHN operators that are generated in the presence of the corresponding mediators. In addition, we define Scenarios 1 and 2, which correspond to consider only right-handed operators, with and without the SM-like contributions, respectively. The set of scenarios that we are going to analyse and the operators involved in each case are: Scenarios 3, 6 and 8 do not generate any left-handed operator, making the "a" and "b" labelling unnecessary. In Scenarios 6, 7a and 7b, where scalar and tensor couplings arise at the NP scale, the renormalization-group running between Λ NP ∼ 1 TeV and the scale m b generates the factor r ≈ 2. Scenarios 3 to 7 have been also studied at Ref. [35]. Within each scenario we will perform a standard χ 2 fit to the data. There are 60 experimental degrees of freedom (d.o.f.), 4 corresponding to R D ( * ) ,F D * L andP D * τ , and 56 to the binned q 2 distributions. Therefore, the number of d.o.f. of our fits is 60 − N WC − 1 = 59 − N WC , where N WC is the number of Wilson coefficients entering in the fit.
All solutions resulting from our fits will present up to three flipped minima with degenerate χ 2 values. The first flipped minimum is obtained by reversing the sign of the LHN Wilson coefficients while keeping the right-handed Wilson coefficients untouched: for X = S, V, T and i = L, R, except for C V LL . The second flipped minimum is obtained reversing only the right-handed coefficients, for X = S, V, T and i = L, R, and the last one flipping both left and right Wilson coefficients, for X = S, V, T and i = L, R, except for C V LL . From now on, we will only discuss the minimum which is closest to the SM scenario.
In the following subsections, we will present the fitted solutions for each considered scenario. Whenever some uncertainties are marked with the symbol † (i.e., C V RR = −0.69 +0.64 † −0.44 ), this indicates that the χ 2 distribution has fallen to another minimum. In these cases, the uncertainty is defined as the range between the central value and the point in which the χ 2 falls to the other minimum. To complete the discussion, it is interesting to see the predicted values of the different observables within each fitted scenario. This information is given in Fig. 9 and in Table 4, where the numerical predictions are marked either with a green tick () if they agree with the experimental value at 1σ or with a red cross () if they do not agree. All minima are in agreement with all experimental observables at the 2σ level.

SM fit
The SM fit, where all the Wilson coefficients are set to zero, i.e. C X AB = 0, gives us the following χ 2 : corresponding to a 69.95% probability (p-value, defined below). The "apparent" good quality of the fit, i.e. The last results suggest an overestimation of the absolute χ 2 value, which is introduced while considering in the fit multiple inputs with large uncertainties as, in our case, the q 2 distributions for the B-meson semileptonic decays. The goodness of a fit is usually characterized through the p-value, defined as where χ 2 (z, n) is the χ 2 probability distribution function with n d.o.f.. Larger p-values correspond to better explanations of the experimental data than lower ones. In order to quantify the quality of our fit, it is convenient to introduce another parameter called Pull that compares any fitted solution with the SM results. This statistical measure is defined as the probability in units of σ corresponding to the difference ∆χ 2 i ≡ χ 2 SM − χ 2 i , assuming that ∆χ 2 i follows a χ 2 distributed function with ∆n i ≡ n SM − n i d.o.f., where the label i refers to the ith scenario. The translation from probability to sigmas is done by associating such probability to the one corresponding to a Pull number of standard deviations in a normal distribution with ∆n i d.o.f., 3 i.e. [70,71] where CDF( In Table 3 we display the Pull SM values of the different fitted minima, together with their corresponding p-values, for all the scenarios analysed. In order to better quantify how favourable are the fitted scenarios with respect to the SM regarding the different observables entering in the fit, we also include their pull for the particular pieces of the χ 2 , splitting it into three contributions: the polarization observablesP D * τ and F D * L , the ratios R D and R D * and the q 2 -distributions of the B → D ( * ) τν decay. In the former we ignore the FF contribution to the χ 2 . As we can see in Table 3, all scenarios exhibit a sizeable improvement with respect to the SM p-value.

Scenario 1: ν R + SM-like
Considering only RHN operators and the SM-like contribution, i.e. C V LL , and imposing an upper bound for B(B c → τν) of 10%, we find two different solutions: a global minimum and a local one with a slightly higher χ 2 , i.e.
Shifting the Wilson coefficients up to 1.2σ, the global minimum becomes compatible with a solution in which the only non-vanishing Wilson coefficients are C V LR and C T RR . As it can be seen in Fig. 2, both C V LR and C T RR help to reproduce the experimental value of R D , R D * andP D * τ . ForF D * L it is a combination of several operators that helps. In the local minimum, the dominant contribution comes from C V RR . As it can be seen in Table 4, both minima saturate the B(B c → τν) ≤ 10% constraint. Thus, relaxing it to be up to a 30%, we find and 3 A probability of (68.3%, 95.5%, 99.7%) equals to (1σ, 2σ, 3σ), respectively.
The value of the χ 2 /d.o.f. slightly improves in this case, whereas the scalar Wilson coefficients are further away from the SM limit. In both cases one can see that most of the Wilson coefficients have large uncertainties. This can be understood from the fact that a large set of variables to fit allow for larger correlations among them, which in turn allows wider ranges for the Wilson coefficients considered. The global and local minima have in fact quite close values of χ 2 /d.o.f., and the χ 2 distribution in the region between them is rather flat. Thus, when evaluating their 1σ variations, one minimum falls often into the other one, as indicated by the † symbols.
This scenario is the most general, in the sense that the preferred C V LL solution without considering RHNs [31] is included in the fit, together with all possible contributions generated as a consequence of having RHNs. No specific NP scenario has been assumed in here.

Scenario 2: ν R
In this scenario we consider solely the contribution to b → c processes coming from the presence of RHNs in the theory. Again, this assumption is very general and model independent, in the sense that no specific types of NP mediators are assumed.
As in the previous scenario, with the constraint B(B c → τν) ≤ 10%, a global and a local minimum are obtained: By shifting all the Wilson coefficients within their 1σ uncertainties, the global minimum is compatible with a solution in which the only non-zero coefficient is C V LR . This coincides with the fit dealing only with the LHN operators where the global minimum was compatible with a global shift of the SM-like operator (i.e. C V LL = 0 ) [31]. In other words, C V LR plays a similar role as the ν L Wilson coefficient modifying the SM contribution. In the local minimum, the main contributions to the observables are coming from C V RR . Since the previous fit saturates the leptonic B c decay bound, we list below the minima obtained after relaxing such constraint to B(B c → τν) ≤ 30%: Similarly to the previous scenario, when relaxing the leptonic decay bound, the χ 2 experiences an improvement and the scalar Wilson coefficients further depart from the SM limit.

Scenario 3: V µ
The mediator V µ ∼ (1, 1, −1) only involves interactions with RHN regarding b → c transitions. Note that we call it V µ instead of the usual nomenclature W µ in order to distinguish it from the SU (2) triplet which does couple to the LHNs. Therefore, this scenario induces exclusively b → cτν R interactions, and particularly the V µ only contributes to the vector Wilson coefficient C V RR .
The global fit gives us the minimum value for this Wilson coefficient together with its χ 2 : Given that in this case our model depends on a single Wilson coefficient, we can study the regions of the parameter space that reproduce the different experimental observables included in the global fit from a fitindependent perspective, as shown in Fig. 3. This figure shows that no region of common overlap can be found at 1σ. This agrees with Fig. 2, which showed that the shift of a single Wilson coefficient with respect to the SM scenario does not modify theF D * L prediction. We also indicate in Fig. 3 the parameter space allowed when relaxing the experimental constraint onF D * L to 2σ and taking B(B c → τν) ≤ 30%. As expected, in that context we find full agreement with the experiment. There is no allowed region forF D * L at 1σ. The lighter orange and red shaded areas correspond to the more relaxed 30% bound on the leptonic B c decay and the 2σ region forF D * L , respectively.

Scenario 4a: Φ
Considering that the mediator Φ ∼ (1, 2, 1/2), with the same quantum numbers as the SM Higgs, is responsible for the NP interactions, and assuming that only right-handed Wilson coefficients appear at the low-energy scale, two different minima with the same χ 2 value, are found. As one can see, they correspond to degenerate solutions, flipping the values of C S LR and C S RR . This can be easily understood by looking at the expressions of B → D and B → D * listed in Eqs. (11) and (20), respectively. These observables depend on the absolute values of the right-handed scalar and pseudoscalar combinations of Wilson coefficients when the vector coefficients are switched off, and therefore remain invariant under the exchange C S LR ↔ C S RR . The same is true for the D * polarization observables that, as shown in Eqs. (34) and (35), are blind to a sign flip of the combination C S RR −C S LR . As Table 4 shows, these minima saturate the B(B c → τν) ≤ 10% bound. Relaxing this constraint to B(B c → τν) ≤ 30%, the minima read where, as expected, the pseudoscalar combination of Wilson coefficients increases its value and the χ 2 slightly improves.
In the left panel of Fig. 4 we show the two-dimensional parameter space where the different observables entering in the fit are satisfied at 1σ. As the figure shows, there is no overlap at this given probability. In this case, not even relaxing the leptonic B c decay upper bound to 30% and the F D * L experimental measurement to 2σ, an overlap in the parameter space is achieved.

Scenario 4b: Φ
The Two Higgs Doublet Models are the simplest examples of UV physics generating this scenario. In addition to RHN operators, a second scalar doublet with the same quantum numbers as the SM one generates LHN Wilson coefficients. The preferred solution of this scenario corresponds to vanishing right-handed Wilson coefficients, which eliminates the degeneracy under C S LR ↔ C S RR . Owing to the interference with the SM-like contribution, an analogous symmetry does not exist for the left-handed coefficients and, therefore, we find in this case a single solution with B(B c → τν) ≤ 10%: With the relaxed limit B(B c → τν) ≤ 30%, the splitting between scalar operators is larger and the χ 2 slightly improves: The right panel of Fig. 4 shows the two dimensional parameter space where the observables entering in the fit are satisfied at 1σ. In this figure, the LHN operators are fixed at their best-fit values. As it can be seen, there is no overlap at this given significance level. The non-existing overlap is also reflected in Table 4 and Fig. 4, where one can see that scalar solutions cannot satisfy R D * , norF D * L . The later is also shown in a very intuitive way in Fig. 2.

Scenario 5a: U 1µ
The presence of the vector leptoquark U 1µ ∼ (3, 1, 2/3) at the high-energy scale will contribute to both left and right-handed operators at the m b scale. This vector leptoquark can be UV-completed in Pati-Salam based unification theories [72][73][74][75] for instance. Considering only the RHN operators, the preferred solution is compatible with a non-zero value of C V RR while C S LR = 0 at 0.4σ, i.e.
Since the scalar coefficient is suppressed, the B(B c → τν) limit is not saturated. Furthermore, all the observables included in the fit agree at 1σ, exceptF D * L which is compatible with the experimental value at 2σ, as illustrated in the left-panel of Fig. 5.

Scenario 5b: U 1µ
Including the contributions to LHN operators, the value of the χ 2 remains almost constant with respect to Scenario 5a, ∆χ 2 = −0.02 for 2 new d.o.f., and the left-handed Wilson coefficients are compatible with zero within 1σ: This indicates that the best solution for a leptoquark with these quantum numbers involves only RHN operators.
The right panel in Fig. 5 shows the small changes on the allowed regions, in comparison with Scenario 5a (left panel). Again, all observables are satisfied at 1σ, except forF D * L .
In this case, there is only one free parameter, since the two relevant coefficients, C T RR and C S RR , are correlated by the Fierz identities. Therefore, one can study the predictions of the fitted observables as a function of only one free parameter in a fit-independent manner, as we show in Fig 6. The region with larger overlap in this figure corresponds to the minimum listed in Eq. (59) and its flipped solution. As in previous scenarios, it is not possible to reproduce the experimental value ofF D * L at 1σ. However, agreement can be find when B(B c → τν) ≤ 30% andF D * L is considered at 2σ. There is no allowed region forF D * L at 1σ. The light orange and red shaded areas correspond to the more relaxed 30% bound on the leptonic B c decay and the 2σ region forF D * L , respectively.

Scenario 7a: S 1
The scalar leptoquark S 1 ∼ (3, 1, 1/3) is considered in this scenario. For Scenario 7a we obtain a solution dominated by a single Wilson coefficient, C V RR , being C T RR compatible with zero within 1σ: The left panel of Fig. 7 shows the regions in the two-dimensional parameter space where the experimental observables can be reproduced at 1σ. Again, at this level of precision, the longitudinal D * polarization cannot be accommodated together with the other measurements, although it is possible to find overlap between all experimental data when the value ofF D * L is taken at 2σ, shown in the figure as a red grid.
For the RHN coefficients, C V RR and C T RR , the χ 2 distribution turns out to be very flat between the two flipped minima, which no longer can be separated. This implies a very broad negative 1σ interval for C V RR , reaching its flipped minimum C V RR = −0.367. As in the case of the vector leptoquark U µ 1 (Scenarios 5a and 5b), the preferred solution for an S 1 leptoquark involves only RHN operators.

Scenario 8:Ṽ µ 2
This is another genuine scenario of RHNs, since it does not generate any b → c transition involving ν L operators. The vector leptoquarkṼ µ 2 ∼ (3, 2, −1/6) only contributes to the Wilson coefficient C S LR . This allows us to study the parameter space preferred by the experiment from a fit-independent point of view. As Fig. 8 shows, there is no overlap among the different experimental constraints at the 1σ level, nor even considering a more relaxed 30% bound for the leptonic decay B(B c → τν) and the experimental value of F D * L at 2σ. Numerically, the fit provides the following minimum:  Table 3 summarizes the fit quality of the results obtained in the different scenarios analysed, quantified through the corresponding χ 2 /d.o.f., the pull with respect to the SM, and the p-value. The resulting predictions in each scenario for the observables included in the fit are also given in Table 4, and compared with their experimental measurements in Fig 9. Several conclusions can be extracted from these results:

Comments on the fit results
• In general, it is difficult to reproduce the experimental value of the longitudinal D * polarization within its 1σ range. From Fig. 9 and Table 4 we can see that the only solutions reproducing all the experimental values (marked with a ) are Scenario 1a with either a 10% ( Min1) or 30% (Min1 and Min2) upper limit on B(B c → τν), and Scenario 4b with a 30%.
• All solutions exhibit pulls between 1.2 and 3.7 with respect to the SM fit, showing a clear preference for NP contributions.
• The largest pull with respect to the SM fit is obtained in Scenario 3, which only contributes to the C V RR coefficient. Note that C V RR plays a similar role than C V LL in the observables involving b → c transitions. Therefore, the preference of the fit for this scenario can be easily understood, since a SM-like modification was the best fit solution in absence of RHN [31].
As Table 4 and Fig. 9 show, Scenarios 4a, 4b and 8 fail badly reproducing the experimental value of R D * .
• Scenarios 4a, 4b, 6, 8 and Scenario 2 Min2, are disfavoured by the q 2 differential distributions of the B → D ( * ) decay with respect to the SM, as the corresponding Pull SM in Table 3 shows.
• Those solutions further away from the SM (larger pulls) present higher p-values, as Table 3 shows.
• In scenarios with several operators, the best fits correspond to solutions where all Wilson coefficients but one are compatible with zero. The non-zero Wilson coefficient is typically C V RR (Scenarios 5a, 5b, 7a and 7b).
• When scenarios with and without LHN operators ("b" and "a" variants, respectively) are compared, the fit indicates a preference for solutions with all left-handed Wilson coefficients compatible with zero within 1σ.
Comparing our results with similar fits previously done in the literature, we can quantify the impact of adding the differential q 2 distributions and considering recently measured observables such asF D * L or P D * τ , together with the update of some experimental measurements. Ref. [35] analysed all mediators that can contribute to the b → cτν R transition, except theṼ µ 2 vector leptoquark, but only included in the fit the values of R D and R D * . The global minimum obtained in Ref. [35] for an extra gauge boson V (Scenario 3) agrees with ours, while the two minima obtained for our Scenario 4a deviate more from the SM solution than ours. The latter is due to the fact that the B c → τν constraint, which has a strong impact on solutions involving scalar Wilson coefficients, was not taken into account in the fit. Indeed, Fig. 2 from Ref. [35] shows that their minima are excluded by this constraint, and this is the reason why in our analysis, this χ 2 is the most unfavorable among all the scenarios considered. For our Scenario 5a, mediated by U µ 1 , two minima are observed in Ref. [35] where the furthest one from the SM solution is ruled out by the constraint B(B c → τν) ≤ 10%. This situation is repeated in the scenario mediated by S 1 , Scenario 7a. Finally, in the case of theR 2 mediator (our Scenario 6), both minima differ slightly from ours since, again, as their Fig. 2 shows, they are excluded by the B c leptonic decay limit; however, taking into account the minimum value of the χ 2 satisfying this constraint, our result is compatible with Ref. [35].  The * symbol indicates that the χ 2 of a given scenario is greater than the SM one.

Predictions
In this section we show the predictions of different observables for the fitted scenarios considered in the previous section. As we will discuss in the following, these results can be used to discriminate between the different scenarios and, in some cases, even distinguish the contribution originated by light RHNs from the SM one.

Predictions of integrated observables
In Table 4 we list the predictions of the different integrated observables considered in the fit, i.e. R D , R D * ,F D * L ,P D * τ and the leptonic branching fraction B(B c → τν), for each of the scenarios considered. Those predictions that are in agreement with the measured values at the 1σ level are marked with a , while a mark indicates disagreement. Only in Scenarios 1 and 4b it is possible to simultaneously satisfy all experimental constraints. The second column shows that the upper bound on the B c leptonic decay is always saturated in Scenarios 1, 2, and 4, which denotes that larger pseudoscalar and axial combinations of the Wilson coefficients would still be preferred.
S7b S8 Figure 9: Predictions for the fitted observables, normalized to their measured values, with their 1σ experimental uncertainties shown as orange bands. For these predictions B(B c → τν) ≤ 10% is taken. The green and red regions indicate the predictions arising from each NP scenario that are in agreement or not with the experimental value, respectively, at the 1σ level. The labels within brackets specify the minimum within a given scenario. The numerical values of these predictions are listed in Table 4.

Predictions of angular coefficients
The three-body differential distribution in B → Dτν and the full four-body angular analysis of B → D * τν → (Dπ)τν provide a multitude of observables that could be experimentally accessible. The presence of neutrinos in the final state makes the measurement troublesome, compared to the case of well-known neutral-current transitions like B → K * µμ. Nevertheless, measuring the distribution of the secondary τ decay, some information on the angular coefficients J i and I i , defined in Eqs. (7) and (14), could be obtained in the near future. As it can be seen from their explicit analytic expressions in Eqs. (9) and (16), these q 2 -dependent functions can be very sensitive to the NP Wilson coefficients present in the theory. In this section, we provide the predictions of such observables in some relevant NP scenarios considered in this work. Fig. 10 shows the predictions for the forward-backward asymmetries A D ( * ) F B defined in Eqs. (12) and (27), the lepton polarization asymmetries of Eqs. (13) and (31) and the longitudinal D * polarization F D * L defined in Eq. (24), as functions of q 2 . For simplicity we have illustrated the four NP scenarios with largest pulls with respect to the SM. Note that Scenario 3, which contains the single Wilson coefficient C V RR , will always give the same predictions as the SM scenario for the forward-backward asymmetries, F D * L (q 2 ) and the angular coefficientsĪ i (q 2 ). Therefore, this scenario is only included in the τ polarization asymmetries. Error bands in these plots correspond only to the uncertainties arising from the fitted Wilson coefficients. These uncertainties have been obtained by minimizing the χ 2 , imposing O i = O i,min + ∆O i,min , and taking the value of the observable O i for which χ 2 = χ 2 min + 1. Other smaller errors such as FF parameters or additional inputs are not taken into account. Therefore the SM predictions, plotted as dotted black lines, do not present any uncertainties.
From these plots, we can see that scenarios with a larger number of Wilson coefficients also have larger uncertainties (Scenario 1, Min 1), as expected because of the wider allowed range of variation of their Wilson coefficients. The forward-backward asymmetry A D F B could be useful to distinguish Scenario 6a from the SM, but the large uncertainties make difficult to discriminate it from other scenarios or to differentiate the SM from Scenarios 1, 6a and 7. A precise measurement of A D * F B would allow to distinguish Scenarios 1 and 6a from the rest of NP scenarios, which partly overlap with the SM prediction. A similar situation occurs for F D * L , where clear differences manifest at low values of q 2 while the different scenarios considered tend to overlap at high q 2 . The τ polarizations P D ( * ) τ are useful to distinguish Scenario 3 from the SM, since these are the only observables that are sensitive to a single shift in C V RR . Moreover, in Scenario 1 P D τ and P D * τ exhibit a quite different dependence on q 2 compared to the other scenarios, which could be exploited to distinguish it at low q 2 values. In the high q 2 region, P D * τ also allows to discriminate Scenario 1 from the other possibilities.
In Fig. 11 we plot the B → D * τν angular coefficients, as functions of q 2 , normalized by the decay width:Ī The CP-odd quantities I 7 , I 8 and I 9 are identically zero in our case, because we have only considered real Wilson coefficients in our fits. It is interesting to notice that despite the large uncertainties Scenario 1, Min 1 can be easily distinguished from the SM predictions and from other minima (for instance looking atĪ 1s orĪ 5 ). However, being able to distinguish other scenarios would be more complicated, unless the current errors on the Wilson coefficients are sizable reduced. There is always an overlap between the SM predictions, Scenario 7a and Scenario 5a. Scenario 6a is close to Scenarios 5a, 7a and the SM predictions, but it is still possible to distinguish it looking at low (Ī 1s ,Ī 5 ) or high (Ī 2s ,Ī 2c ,Ī 3 andĪ 4 ) q 2 values.

Conclusions
Using an EFT approach, we have explored the impact of various NP operators on the recently observed anomalies in b → cτν transitions. In particular, the focus of this work has been to identify the role of NP operators which can arise due to the presence of RHN in the theory. This has been achieved through a globalfit analysis of all available b → cτν data until date: R D ( * ) ,P D * τ ,F D * L and the q 2 differential distributions of B → D ( * ) . Previous analyses only studied the integrated rates and did not include the polarization information (P D * τ ,F D * L ) and the q 2 distributions measured by the BaBar and Belle collaborations, which play an important role in discarding many proposed NP explanations.
We have also studied the differential B → Dτν decay distribution and have derived the full fourbody angular distribution of the decayB → D * (→ D π)τν, for the most general dimension-six effective Hamiltonian, which includes (axial)vector, (pseudo)scalar and tensor operators for both the left-and righthanded leptonic currents. The rich dynamical information embodied in the coefficients of these angular distributions could be, in principle, experimentally accessed. From these distributions, we have constructed different observables and have analysed their predicted values within the NP scenarios emerging from our fits. In the next few paragraphs, we briefly summarize the key findings of our analysis.
NP contributions have been assumed to be present only in operators involving charged leptons of the third generation, which is well justified since potential NP effects in b → c ν transitions ( = e, µ) are known to be negligible [18]. The NP couplings have been also assumed to be real, due to the absence of any evidence of CP violation in these channels. After investigating the separate impact of individual Wilson coefficients, we have performed multi-dimensional fits to the data within eleven different scenarios. The first and the second case include all five RHN operators with and without a SM-like NP contribution, respectively, whereas the remaining scenarios correspond to 'simplified models' obtained by integrating a single mediator above the EW scale: namely, a scalar boson Φ, a vector boson V µ , two scalar leptoquarks S 1 andR 2 , and two vector leptoquarks U µ 1 andṼ µ 2 . In those cases where the tree-level exchange of a mediator generates both ν L and ν R operators, we have further analysed two model variants with and without the ν L contributions.
Among all scenarios analysed, the vector boson V µ (Scenario 3) seems to be the preferred option, in terms of the pulls from the SM hypothesis, as shown in Table 3. The next two possibilities are the scalar leptoquark S 1 (Scenario 7) and the vector leptoquark U µ 1 (Scenario 5), switching on the RHN couplings only, which can also provide good agreement to the data. However, it is important to note that none of these three possibilities can generate values of the longitudinal D * polarization within its current 1σ experimental range; they can only reach agreement with theF D * L measurement at the 2σ level. Interestingly, theF D * L data can only be explained at 1σ in very few cases, namely, with all RHN operators plus the SM-like contribution (Scenario 1), or with a scalar boson Φ, switching on both ν L and ν R operators (Scenario 4b) and with a relaxed upper limit of 30% on B(B c → τν). However, these scenarios are not the best choices in explaining the R D ( * ) measurements in terms of pull, as reflected in Table 3. Nevertheless, they do reduce the R D ( * ) deviation significantly, and bear very important information about simultaneous agreement of all observables considered in this work. Due to the large uncertainty of the currentP D * τ measurement, all scenarios are compatible (within ±1σ) with it. The R D measurement is also easily accommodated in all the NP scenarios that we have analysed.
Measurements of additional observables such as polarizations and angular distributions could help to to disentangle the dynamical origin of the current anomalies. In particular, we have displayed the information contained in the three-body and four-body angular distributions of B → Dτν and B → D * (→ Dπ)τν, respectively, and their sensitivity to the different NP scenarios analysed. The experimental measurement of these distributions is of course very challenging because of the presence of undetected neutrinos, and one would need to further analyse the decay products of the tau in order to recover the accessible information.

Acknowledgements
This work has been supported in part by the Spanish Government

A Form factors
The hadronic matrix elements can be parametrized by showing explicitly their Lorentz structure as [31,76] for the process involving B → D, and the B → D * hadronic matrix elements as The FFs F 0 (q 2 ), F 1 (q 2 ) and F T (q 2 ) appearing in the B → D matrix elements are defined as while the B → D * helicity amplitudes involve the following FFs for vector, axial and pseudoscalar currents, and for the tensor matrix elements, The reduced functionsĥ i (q 2 ) = h i (q 2 )/ξ(q 2 ) take the form [15] h + = 1 +α s C V 1 + ω + 1 2 for B → D, and for B → D * . The explicit expressions of the ω(q 2 )-dependent factorsL 1...6 and the O(α s ) corrections C i can be found in Ref. [15]. Note that corrections of order Λ 2 QCD /m 2 c are included via the subleading Isgur-Wise functions l 1,2 (ω). The detailed parametrization of the different FFs can be found in Ref. [15,18].

B Helicity Amplitudes
We compile here the whole set of relevant helicity amplitudes, following the standard formalism for semileptonic B decays [47][48][49].

B.1 Leptonic amplitudes
The leptonic helicity amplitudes are defined as, (71) where λ τ denotes the sign of the tau lepton helicity and µ (λ) are the polarization vectors of the intermediate virtual boson (λ = t, 0, ±) in its rest frame, defined in Appendix A of [47]. Notice that the helicity of the neutrino is explicitly given in the above equation with the symbols L (λ ν = −1/2) and R (λ ν = +1/2) for left-handed and right-handed neutrinos.

B.2 Hadronic amplitudes
The hadronic helicity amplitudes of B → M τν τ (M = D, D * ) transitions, H λ M i,λ , are defined through the matrix elements [76] and for B → D * :

B.3 Total amplitude
Amplitudes corresponding to left and right-handed neutrinos do not interfere since they correspond to different final states. Using the completeness relation of the polarization vectors µ (λ), λ δ λ * µ (λ) ν (λ) = g µν with δ 0 = δ ± = −δ t = −1 , the decay amplitudes for the transitions B → D ( * ) (λ M ) τ (λ τ ) ν X with X = L, R can be written as where λ τ = ± 1 2 denotes the τ helicity in the rest frame of the τ ν pair. For the D * final state, λ M = 0, ±1 is the D * helicity in the B rest frame, while λ M = s labels the corresponding D pseudoscalar meson. Note that the terms proportional to C T LR or C T RL vanish identically. The helicity amplitudes for the transition B → Dτ ν are defined as where λ τ = ± 1 2 denotes the τ helicity in the rest frame of the τν pair. The decay B → Dτν is characterized by four reduced amplitudes M s,λτ L,R , which are given by: The q 2 -dependent functionsÃ X λ (λ = 0, t, S, T ; X = L, R) have been defined in Eq. (10). Analogously, we express the B → D * τν transition in terms of twelve reduced helicity amplitudes while the corresponding amplitudes for RHNs are given by: The transversity amplitudes A L,R λ are listed in Eq. (17).

C Propagation and decay of D *
To compute the full four-body decay amplitude B → D * τν → (Dπ) τν, we need to describe the propagation and the decay of the vector boson D * to the Dπ final state. The D * → Dπ amplitude can be parametrized in the form M with an effective coupling g D * Dπ that can be determined from the total decay width, where C = 1, 1 2 for a final π ± , π 0 , respectively. The dependence of the effective amplitude (86) on the momentum and polarization vectors fixes the angular structure of the three possible helicity amplitudes: with | p D | = λ 1/2 (m 2 D * , m 2 D , m 2 π )/(2m D * ) being the three-momentum of the D meson in the D * rest frame. The propagation of the D * can be described through a Breit-Wigner function. Since the decay width of the D * is much smaller than its mass, we can use the narrow-width approximation, and write the decay probability of the process B → (Dπ) τν in the form Notice that the dependence on g D * Dπ cancels out from this expression. The interferences among the unobservable helicity amplitudes of the intermediate D * meson generate the different dependences on θ D , appearing in the four-body angular distribution listed in Eq. (14).