$b \to c \tau \nu_{\tau}$ Decays: A Catalogue to Compare, Constrain, and Correlate New Physics Effects

In this article, we have predicted the standard model (SM) values of the asymmetric and angular observables in $B\to D^{(\ast)}\tau\nu_{\tau}$ decays, using the results of the new up-to-date analysis in $B\to D^{(*)}\ell\nu_{\ell}$. We have also revisited the SM prediction of the inclusive ratio $\mathcal{R}_{X_c}$, and have given its values in different schemes of the charm quark mass. This is the first analysis which includes all the known corrections in the SM. In addition, we have analysed the $b\to c\tau\nu_\tau$ decay modes in a model-independent framework of effective field theory beyond the standard model. Considering all the possible combinations of the effective operators in $b \to c \tau\nu_{\tau}$ decays and using the Akaike Information Criterion, we find out the scenarios which can best explain the available data on these channels. In the selected scenarios, best-fit values and correlations of the new parameters are extracted. Using these results, predictions are made on various observables in the exclusive and inclusive semitaunic $b \to c $ decays. The graphical correlations between these observables are shown, which are found to be useful in discriminating various new physics scenarios.


Introduction
The semileptonic B → D ( * ) ν ( = µ or e) decays play important role in the extraction of |V cb | as well as the extraction of the form factors associated with the b → c decays. These form factors are the major inputs used in the predictions of R D ( * ) = B(B → D ( * ) τ ν τ )/B(B → D ( * ) ν ), where = µ or e; for recent updates see [1][2][3][4][5][6][7]. There are additional form factors which can not be extracted directly from the experimental data and one thus needs to rely on the heavy quark effective theory (HQET) inputs. On the other hand, the lattice simulations can predict these form factors at zero and non-zero recoils [8,9]. The standard model (SM) predictions of R D ( * ) which are used in the literature are given by [1,10] R D = 0.299 ± 0.003, R D * = 0.252 ± 0.003. (1) The prediction of R D includes the up-to-date lattice inputs, whereas the prediction of R D * relies heavily on HQET inputs. Also, the current lattice results suggest that the HQET values of form factors at zero recoil are not in complete agreement with those from lattice [1,3]. This discrepancy could be due to the missing higher order (α 2 s and α s Λ QCD /m b ) corrections in the HQET relations of the form factors. The inclusion of lattice inputs increases the value of R D * [2,5,6]. Moreover, the data allow the Send offprint requests to: unknown corrections in the ratios of the HQET form factors to be as large as 20% [6]. With all these inputs and using the Caprini-Lellouch-Neubert (CLN) [11] parametrisation of the form factors, the SM prediction is given by R D * = 0.259 ± 0.006 [6]. On the experimental side, the current world averages are given by [12] R D = 0.407 ± 0.046, R D * = 0.304 ± 0.015 . (2) We note that the deviations in both the observables are a little less then 2.6σ. Though these deviations can be explained by a variety of new physics models, we will follow a model independent analysis, like the one done in ref. [13]. The other model independent analyses, which are relatively new, can be seen in references. [14][15][16][17][18][19][20][21][22]. There are a few other observables, which could be constructed from the B → D ( * ) τ ν τ decays and which are potentially sensitive to the new physics (NP) beyond the SM. For an update, please see [23] and references therein. Among these, the τ polarisation asymmetry has been measured by Belle [24], and it is consistent with the SM predictions (see section 2 for detail).
Recently, LHCb has published their result on another b → cτ ν τ decay mode, where the observable and its measured value are given by [25] = µ or e. Though the uncertainties are large, this measurement is 2σ above the corresponding SM prediction, which lies in between 0.25 and 0.29 [26][27][28][29]. We note that in both R D ( * ) and R J/ψ , the measured values are above the SM predictions and therefore, the NP should contribute constructively to both the decay modes, in order to explain the observed discrepancies. There are a few other b → cτ ν τ decay modes like B c → τ ν τ , Λ b → Λ c τ ν τ and the inclusive decay B → X c τ ν τ which are potentially sensitive to the new interactions and the NP affecting R D ( * ) and R J/ψ should also have an impact on these decay modes.
Therefore, the correlation studies of the various observables associated with these decay modes will be an important probe for an indirect detection of NP. On the other hand, the precise measurements of some of these observables will be useful to constrain the new physics parameters associated with a model. This motivates one to predict the values of all the relevant observables for some specific model, which can then be further checked for consistency with the future measurements.
In this article, we have done a model independent analysis of the NP affecting the b → cτ ν τ decay modes. The operator basis is exactly the same as that given in our earlier works [13,23], which consists of scalar (S), vector (V), and tensor (T) type of operators. We did not consider the scenarios with right handed neutrinos. We have considered all possible combinations of these operators and categorised them as independent models. There are several models capable of describing the observed data and one is thus confronted with the problem of model selection. We use the Akaike information criterion (AIC) to find out the best possible model(s) for the existing data. A model selection criterion is a formula that allows one to compare models; for details, see [30]. Alternative related approaches to the model selection are the bootstrap method and crossvalidation. Cross-validation works poorly with small sample sizes, as it is in our case, and parametric bootstrap variants of AIC have recently been proposed [31], which we have not used in the present analysis. Using the AIC, we have first selected the models best suited for explaining the existing data. Then, with the best fitted values of the model parameters, we have predicted the values of various observables associated with the above mentioned decay modes. We have studied the correlations amongst the observables in detail as well.
2 New Physics and the observables 2

.1 New operators
In this subsection, we will discuss the complete operator basis in b → cτ ν τ decays. The most general effective Hamiltonian describing these transitions is given by where C W (W = V 1 , V 2 , S 1 , S 2 , T ) are the Wilson coefficients corresponding to the following four-fermi operators: Here, we have considered only the left handed neutrinos.

Observables in b → cτ ν τ decays
We will define various observables used in our analysis in this subsection.

B → D ( * ) τ ν τ
Following ref. [13] and references therein, we can write the differential decay rates for B → D ( * ) τ ν τ with the Hamiltonian in eq. 4. The q 2 -distribution of the decay rate of the decays B → D ( * ) ν are obtained by setting C W = 0 and m τ = m . Here, we are assuming that new effects are present only in B → D ( * ) τ ν τ and the B → D ( * ) ν channels for = µ and e are free from any NP effects. H s V,Y (q 2 ) and H V,Y (q 2 ) are the helicity amplitudes forB → D and B → D * transitions respectively (with Y = ±, 0 and t). These amplitudes can be expressed in terms of form factors in B → D ( * ) transitions. The details of the form factors and their parametrisations can be looked up in [6] and the references in there.
In the present work, we have followed the CLN [11] parametrisation of the B → D ( * ) form factors and have used both the fitted and predicted values of these parameters obtained in [6]. In terms of the differential distributions, the ratios R D ( * ) are defined as with q 2 max = (m B − m D ( * ) ) 2 , and = e or µ. Along with these ratios, there are a number of other observables, that can be constructed in these channels, which are sensitive to NP. Most of them are not yet measured experimentally. These are: τ -polarization is defined by studying further τ decays: where , λ τ is the τ helicity, and q 2 max = (m B − m D ( * ) ) 2 .
-D * longitudinal polarization can be extracted from the angular distribution in D * → Dπ decays: where -If we write the double-differential decay distribution as where θ is the angle between the three-momenta of τ andB in the τν rest frame, then b ( * ) θ (q 2 ) determines the lepton forward-backward asymmetry in the following way: As mentioned earlier, in addition to these observables there are several other channels that will be affected by the same set of NP operators. We have used some of those most relevant observables in our analysis, either as fit inputs or as constraints and/or for prediction.

B c → J/ψ ν
Ratios similar to those defined in eq. 6 can be defined for the decay channelB c → J/ψ ν by replacing the respective mesons. For various form factors in B c → J/ψ decays, see [32]. Given the unavailability of a precise calculation of B c → J/ψ form factors till date, we have the option to choose from a collection of available parametrisations [27][28][29][33][34][35][36]. Choosing different parametrisations results in varying the central value of R J/ψ within the range 0.25 -0.29, which is considered theoretical range in recent experimental analyses. Taking the uncertainties from different parametrisations into consideration, we see that the allowed theoretical range of R J/ψ is actually larger than that. We consider two parametrisations residing at two far ends of this range, namely perturbative QCD (PQCD [29]), and light-front covariant quark model (LFCQ [36]) in this work. A preliminary result on the form factor A 1 (q 2 max ) (this is the only form factor contributing to the decay at zero recoil) is available from the HPQCD collaboration [37]. This result is consistent with the parametrisations used in this draft.
The q 2 distribution for the decay process (Λ b → Λ c τ − ν τ ) can be written as [46] dΓ where represent the contributions from the vector and axial vector currents respectively. Their origin could be either the SM or any NP model. A SP 3 and A T 4 represent the contributions from the scalar-pseudoscalar and tensor currents, which will appear only in the NP models.
are the interference terms which will have contributions from various operators in the SM, as well as an NP model. These are functions of combinations of the helicity amplitudes H λ Λc ,λw , which in turn can be expressed in terms of form factors and NP couplings. Several instances, where these form factors have been studied using sum rules and quark models, can be found in the literature [38,[47][48][49][50][51][52][53][54][55][56][57]. For our purpose, helicity form factors have been calculated using the formula from lattice QCD in the relativistic heavy quark limit [39].
Similar to the ratios defined earlier, two observables can be defined here, motivated by the lepton flavor university violation elsewhere: Along with these ratios, we have also considered the forward-backward asymmetry in Λ b → Λ c τ − ν τ , defined as where Γ (1) = dΓ d cos θτ and θ τ is the angle between the momenta of the τ lepton and Λ c baryon in the dilepton rest frame.

B → X c τν τ
Similar to the ratios R D ( * ) , we can define the following ratio for the inclusive decay B → X c τ ν τ : with = µ, e. The decay B → X c ν is well studied in the literature; for a comprehensive update see [58] and references therein. In the present work, the detailed mathematical expression of the decay width of B → X c ν and all    The correlations between various non-perturbative parameters and the masses. These were all obtained in the analysis of inclusive B → X c ν decays in [40].  Table 4: Present experimental status of the observables used in this analysis. First uncertainty is statistical and the second one is systematic. * This correlation is between R(D * ) and Pτ (D * ). Stat. corr. = 0.29 and syst. corr. = 0.55. † This uncertainty originates from the uncertainties on B(B 0 → D * − π + π − π + ) and B(B 0 → D * − µ + νµ). other relevant inputs are taken from [40]. The simplified expression for the decay width of the inclusive semitaunic decay of B meson in SM are given in [59]: Here, the terms involving C 0 , and C (2) 0 represent the contributions from the leading order(LO), next-toleading order (NLO) [60], and next-to-next-to-leading order (NNLO) [61] corrections in α s respectively, whereas C µ 2 π , C µ 2 G , and C ρ 3 D , C ρ 3 LS are the contributions at order 1/m b 2 [62] and 1/m b 3 [59], respectively. These coefficients depend on the quark and lepton masses and Γ 0 , defined as The parameters like µ 2 π , µ 2 G , ρ 3 D , ρ 3 LS are the matrix elements of the operators of dimension five and six, respectively, which are non-perturbative in nature. We have also included the well known electroweak correction A ew (= 1.014). As mentioned earlier, the values of the various nonperturbative parameters and the masses (eq. 16) along with their correlations are taken from tables II and III of ref. [40]. In our analysis, the b quark mass is defined in the kinetic scheme, while the c quark mass has been defined in both the kinetic (m Kin c = 1.091 (20) GeV [63]) and the M S scheme (m c (3GeV) = 0.9843(56) GeV [64]). Relations of the pole masses (eq. 16) with the kinetic and M S masses are taken from [65] and [66], respectively. We have also considered α s = 0.22 ± 0.018.
We have given the predictions for the ratio R Xc instead of Br(B → X c τ ν τ ). As can be seen from the above expressions, this ratio is relatively clean, since the errors due to |V cb | and the mass of the b-quark cancel in the ratio. Our predictions for R Xc in the SM are given in table 2. These predictions differ from each other due to the difference in the mass of the charm quark in two different schemes. We note that the central values of the two predictions change by ≈ 2% due to scheme dependence, albeit being consistent within 1σ uncertainties. Also, we have checked our prediction for the 1S scheme masses of the b and c quark, and we agree with that given in ref. [?] which is also different from the predictions given in table 2 (NLO and 1 ). These results are clearly scheme dependent.
In the case of m kin c (1 GeV), the correlation matrix for the non-perturbative parameters and the masses are given in table 3, which are obtained from the analysis of [40]. This is the first analysis which includes all the known corrections in the prediction of R Xc . In ref. [59], the analysis has been done with similar set of inputs, without considering the NNLO corrections. We have checked that our result agrees with them, within the error bar, at the same level of accuracy. The inputs for the analysis with m c (3GeV) are taken from table II of ref. [40]. In this scheme, the predictions have larger uncertainties compared to those in the kinetic scheme. This is due to the difference in the correlation matrix of parameters given in table 2.
To calculate the effects of new physics in the inclusive decay B → X c τν τ , we decompose the decay width as Here, the first piece is arising solely from SM, while the second and third terms are the contributions from NP with different powers of the new couplings. The expressions of Γ N P (1) and Γ N P (2) are taken from [67]. Some other recent works, discussing NP effects in the inclusive mode, are ref.s [68,69].
In terms of the general hamiltonian defined in eq. 4, the branching fraction of B c → τ ν τ can be expressed as [70], where f Bc = 0.434(15)GeV and τ Bc = 0.507 (9)

Experimental status
All the experimental results used in the analysis of NP is tabulated in table 4. There have been quite a few measurements of the ratios R D ( * ) in recent years. Apart from the most recent ones measuring R D * , they are consistent with a sizeable deviation from the SM. The experimental result most deviated from the SM predictions is still the first one reported by BABAR. Though it is apparent from the recent measurements that R D * values are coming down towards the SM, it is still too early to consider it as a trend for two reasons: (a) The experimental uncertainties are still quite large. (b) The actual deviation depends heavily on the correlation between R D & R D * , and any analysis bears the risk of being inconclusive without the simultaneous measurement of both of them. As an example, one can check the Belle (2015) result [42], where the R D * is consistent with the SM result within 1σ, but the combined result is at tension with the SM due to R D and its correlation with R D * .
The first, although quite imprecise, measurement of τ polarization asymmetry is done by Belle in 2015 [42]. Though essentially it is an upper limit, we have included this measurement as a data point in our analysis. Table 4 also contains the recent measurement of R J/ψ by LHCb [25]. Not only is this result in tension with the theoretical predictions, the central measured-value is almost double of that predicted by SM. As the experimental uncertainty is large, they are still consistent with 90% C.L. range. LHCb has used a z-expansion parametrization [71] for the shared form factors for the signal and normalization modes and has determined them directly from the data. As is evident from the theoretical results for R J/ψ (table 1), the PQCD result is a little closer to the LHCb result. As has been pointed out in [72], and later also corroborated in [73], if uncertainty decreases but the central value remains approximately the same in future experiments, NP effects which explain the increase in R D ( * ) , will be unable to explain the measured value of R J/ψ . This should, in essence, result in a worse fit while this value is considered. The decay B c → τ ν, despite being out of the experimental reach for now [74], can be used as an effective constraint on any NP effects that could potentially explain the R D ( * ) and R J/ψ excesses. A conservative upper limit quoted for B(B c → τ ν), even after adding NP effects, is 30% [14]. A stronger upper bound of 10% is obtained from LEP data taken at Z-peak [18] with a prospect of an even tighter bound from the full L3 data [75]. In our analysis, we have used these two constraints.

Numerical Optimization
As mentioned earlier, the goal of this paper is to perform a model independent multi-scenario analysis with the experimentally available results on the charged current anomalies, in conjunction with other relevant results, to obtain a data-based selection of a 'best' scenario and ranking and weighting of the remaining scenarios from a predefined set. If we consider the NP Wilson coefficients occurring in eq. 4 to be complex, all possible combinations of the real and imaginary parts of the coefficients (10 parameters in total) should constitute such a predefined set, from which we can choose different scenarios. Scenarios containing only imaginary Wilson coefficients are neglected.
For each such scenario k, we define a χ 2 statistic, which is a function of the real and/or imaginary parts of the Wilson coefficients (C k W ) associated with the scenario in question, and is defined as: Here, O th p (C k W ) are given by eqns. 6, 7 and sec. 2.2.2 as applicable and O exp p is the central value of the p th experimental result. Statistical (systematic) covariance matrices V stat(syst) , are constructed by taking separate correlations, wherever available. The nuisance parameters (Table 9) occurring in the theoretical expressions are tuned in to the fit using a term where I p k and v p k are the k th input parameter and its respective value. For each scenario, we perform two sets of fits. First, we use different combinations of the experimental results of R D ( * ) (and P τ (D * )). For the second set, we redo the fits including R J/ψ . As the form factor parametrization, as well as the single experimental result for R J/ψ are quite imprecise, instead of defining a χ 2 nuis. (R J/ψ ), we add the SM uncertainty of R J/ψ in quadrature to the experimental one. Following the discussion in sec. 2.2.2, we do two sets of fits in this stage, with two different sets of form factor parametrization for B c → J/Ψ , namely LFCQ and PQCD.
After each fit, we determine the quality of it in terms of the p-value obtained corresponding to the χ 2 min values and the degrees of freedom (DoF) for that fit. We also double check the quality of the fit and existence of outliers in the fitted dataset by constructing a 'Pull' for each data-point and checking the normality (i.e. the probability that it is consistent with a Gaussian of µ = 0 and σ = 1) of their distribution. Due to the small number of data-points in this analysis, no readily available normality test can perform with certainty and it is necessary to scrutinize each individual pull distribution. Still, we perform a variant of the "Shapiro-Wik" normality test as an extra criteria for elimination of scenarios. In other words, we drop the fits which have a pull distribution with probability to be a normal distribution ≤ 5%. Finally, we add the constraints according to sec. 2.2.5 and 2.3 to our analysis and obtain the allowed parameter space. Next, we perform a model-selection procedure on the remaining set of viable scenarios for each data-set. In the following sub-section, we elaborate the method used to do the multi-model selection procedure.

Model Selection Criteria
One measure of the degree of structure in a model, inferred from the data, is its dimension, i.e. the number of parameters in it. In general, bias decreases and variance increases as the dimension of the model increases. Requirement of the optimum dimension of the model is called the 'concept of parsimony' [76], which is essentially a bias versus variance trade-off in statistical terms. All model selection methods, to some extent, depend on the principle of parsimony [77].
The most generally applicable and reliable method for model comparison is 'cross-validation', which, in addition to testing the predictive power of the model, minimizes the bias and variance together by minimizing the meansquared-error (MSE). The problem of applying cross validation to the present analysis is that its applicability to very small sample sizes (as is the case here) is questionable [78,79]. Thus, to the goal of model selection, we have used information-theoretic approaches, especially the second order variant of Akaike Information Criterion (AIC c ) [80] in the present work. This criterion and other competing criteria have previously been applied in one work of ours analyzing the charged current anomalies of b-decay [13]. In that analysis, we had worked with the Table 5: Scenarios selected after passing the normality check and the criterion ∆AICc ≤ 4, for all data available (with or without R J/Ψ ). First and second columns of each dataset represent the reduced χ 2 and corresponding p-value. Third, fourth and last columns represent the independent fit parameters, Akaike weights, and whether or not the fit results satisfy the constraint B(Bc → τ ντ ) ≤ 30% respectively. ' !' means that only some of the multiple minima satisfy this limit for the scenario in question. The scenario where CV 1 is complex (i.e. both Im(CV 1 ) and Re(CV 1 ) present), we consider this as a single parameter case for model selection   binned data on the differential decay distribution of these channels. Given the notation for full reality or truth is f and an approximating model in terms of probability distribution is g, we can create a χ 2 function in terms of the parameters of g and empirical results, following sec. 3.1. For each model g i in a set with R competing scenarios, we can define an AIC c in terms of the χ 2 min in the parameter space which is equivalent to the maximum point of the empirical log-likelihood function: where n is the number of data points and K is the number of estimable parameters 1 in g i . As a rule of thumb, use of AIC c is preferred in literature when n/K < 40. In 1 There is a subtle point here about the number of 'estimable parameters'. As an example, there are two competing scenarios in the present analysis: (a) with non-zero Re(CV 1 ) and (b) with non-zero Re(CV 1 ) and Im(CV 1 ). Now, the two parameters of application, the model with the smallest value of AIC c is estimated to be the 'closest' to the unknown reality generating the data, among the considered models. Whereas all AIC i c are on a relative scale and are strongly dependent on sample size, simple differences of them (∆ AIC i = AIC i c − AIC min c ) estimate the relative expected information loss case (b) always appear together in an identical manner in the expressions of our observables, essentially making the number of 'estimable parameters' = 1. Thus, in the first column of table 5, cases 1-4 has the same DOF = 8. We follow this throughout the analysis. Another way of finding the 'number of estimable parameters' is to calculate the p-value of the fit from toy Monte-Carlo (MC) method. This value, in conjunction to the approximation that the fit-statistic follows a χ 2 distribution, can give us the number of degrees of freedom, and thus the number of estimable parameters (as we have also checked). As we need ∆AICc instead of the absolute value of the AICc in our analysis, the naive way of estimation of the number of parameters (except the case CV 1 , as explained above and which we treat as a special case) works just fine.   Table 9: Nuisance inputs to create χ 2 nuis. defined in eq. 21. These are obtained from the analysis in ref. [6].

Parameters
Value Correlation

Results
Following the methodology described in the previous subsection, we have taken several combinations of the available data and have performed the analysis in the following stages for each dataset.

Model selection
We have created the χ 2 k statistic for the k th scenario containing real and imaginary parts of Wilson coefficients C k W , and repeated that for all k (let us reiterate here that scenarios containing all imaginary C W s are neglected). We have taken scenarios with as many as 4 individual components of Re(or Im)(C W ). Then we have minimized each of those over the corresponding parameter space (with the form factor parameters as nuisance parameters). After checking normality for each fit and dropping scenarios with ≤ 5% significance, we have arranged the remaining scenarios in ascending order of AIC c and have kept only those with ∆AIC c ≤ 4. These are, essentially, the best scenarios to explain the data in that specific dataset under the present experimental constraints.
Tables 5, 6, 7 and 8 contain the listed scenarios of the data sets which are obtained from table 4. Each table essentially contains three variations of similar datasets: the first one is data without R J/Ψ and the rest two, with it. The reason for treating R J/Ψ separately is the apparent tension of the measured central value with that of the SM one, as explained in section 2.3. Moreover, the theoretical values of R J/Ψ are heavily dependent on form factor parametrisation and differ considerably over different choices of it, as explained in section 2.2.2. Thus without showing bias to a particular type of parametrisation, we treat LFCQ and PQCD separately in second and third datasets of each table respectively. As PQCD predicts relatively higher values for R J/Ψ , fits are generally better for these sets than those corresponding to the LFCQ ones, as can be checked from the p-values listed in the second column of these datasets. The measured value of P τ (D * ) has large error. Therefore,we have dropped τ -polarization asymmetry from the list of inputs in one of the fits (table 6). The measured values of R D ( * ) by BABAR are relatively old, and are largely deviated from the respective SM predictions. Therefore, in order to check the impact of BABAR data on our model selections, in one of the fit scenario we have dropped the BABAR data (table 7). We have also done the analysis with the measured R D * alone, which will help us to figure out the sensitivity of this observable towards a particular type of NP scenario.
We notice that for all datasets, the maximum number of independent fit parameters (for listed scenarios) is 2. As is explained in the previous sections, this is natural, because AIC c penalises the increased variance associated with increase in number of independent parameters. The fourth column in these tables, for each dataset, lists the w AICc for each scenario, which estimates the relative likelihood for that scenario (among the given set of scenarios, the number of which is ≈ 90 ∼ 95 for our analysis) to  explain the data. As can be seen, the first few scenarios take up a large chunk of the total likelihood and it is evident that all unselected scenarios together constitute a very small fraction of it. This is another way of understanding why the listed scenarios are the best ones suited to describe the given dataset.
Once we have listed the best scenarios, we scrutinise the allowed parameter space for each of them. As all the finally selected models have a most of two parameters (other than the nuisance parameters), we can plot the marginal confidence levels in the fit-parameter-space with the help of the defined χ 2 function. We have prepared plots for each of these scenarios either in terms of the goodness-offit contours (for two parameter scenarios), or by directly plotting the χ 2 with respect to the parameter (for single parameter scenarios). The contour plots are prepared with constant χ 2 contours equivalent to 1σ and 3σ that correspond to confidence levels of 68.27% and 99.73% respectively. For two parameters, the constant χ 2 values are = χ 2 min +∆χ 2 , where χ 2 min is the minimum value of χ 2 obtained after minimization over the parameter space, and ∆χ 2 = 2.296 and 11.83 for 1σ and 3σ respectively. Sim- Table 10: Best fit results and correlations of the scenarios listed in table 5. We omit the scenarios disallowed by the constraint B(Bc → τ ντ ) ≤ 30%. For the cases where only some of the minima are allowed by the constraint, we quote only those of them allowed by the constraint and closest to the SM point at the same time. Elaboration is in section 4.2. For scenarios where the best fit, instead of being an isolated point, is actually a contour in the parameter-space, we ask the reader to check the corresponding plot. Figure 1c is the plot for scenario 3 in the last dataset of this  (15) ilarly, for single parameter case, 1σ and 3σ intervals are shown with ∆χ 2 = 1 and 9 respectively. As a representative case, we have shown the plots for all scenarios listed in the last dataset of table 5 ('All data', where the theoretical value of R J/Ψ is calculated in PQCD) in figures 1 and 2.
We then use the limits on B(B c → τ ν) mentioned in section 2.3 as our conservative constraint on each scenario. In the aforementioned tables (5, 6, 7 and 8), the last column of each dataset indicates whether the corresponding scenario passes the constraint of B(B c → τ ν) ≤ 30% (' ') or not ('× × ×'). For many scenarios, there are multiple bestfit points (or at least, multiple 68% confidence regions). In some cases, some of these multiple minima are ruled out from the B(B c → τ ν) ≤ 30% constraints, while the rest are allowed. These scenarios are marked as ' !' in the said column. By observing the nature of the confidence levels in the parameter space, we can pick out the minima which are allowed.
In case of the plots depicting the parameter space, we show two limits: regions disallowed by B(B c → τ ν) ≤ 30% (gray shaded region) and B(B c → τ ν) ≤ 10% (diagonally hatched region) respectively. The first one is the conservative limit from B c decay width and the second aggressive one is motivated from the studies quoted in section 2.3 (though we are calling it aggressive, it is a perfectly reasonable upper bound, given the ∼ 2% SM prediction as shown in table 1). We know that even if present, any NP effect is going to be small. Thus when multiple minima are allowed by these constraints, we quote the results for those which are closest to the origin (corresponding to SM) in further analysis 2 . One more thing to note here is that scenarios with either real or imaginary part of C T will not have any constraint on the corresponding axis in the parameter space, as B(B c → τ ν) is unaffected by tensor interactions.
Our results of the model selection with all the available data are shown in table 5. We note that the best possible scenarios are the models with either of the operator O T or O V1 with real C W . The explanations with the operator O S1 are disfavoured by the bound from B(B c → τ ν τ ). These one-operator scenarios are allowed even if we choose ∆AIC c ≤ 2. The same criterion does not allow the twooperator scenarios. However, the constrain ∆AIC c ≤ 4 allows a few of the two-operator scenarios. In all these less preferred two-operator scenarios, the operators like O V2 , O S1 , O S2 appears in combination with either of O T and O V1 . Also, we note that the operator O V2 with complex C W is favoured by the data. Re(CV 1 ) 0.099 ( Re(CS 1 ) 0.178 (67)   Re (   As can be seen in table 6, our conclusions on the selected models will not change if we drop the experimental input on P τ (D * ) from our fit. Also, if we drop the BABAR data on R D ( * ) from the list of the inputs for the fit, the best preferred scenarios are still the one with the operator O T or O V1 with real C W s. However, there are other one-operator scenarios with O V2 or O S1 with real C W s which are then allowed by the criterion ∆AIC c ≤ 4; see table 7 for detail. Here too, there are a few two-operator scenarios which successfully pass the above mentioned criterion. Finally, in order to understand the impact of the R D * , we have also done an analysis considering only the data on R D * (table 8). It shows that all of O V1 , O T , O V2 , O S1 and O S2 are allowed by the criterion ∆AIC c ≤ 4. However, the constraints from B(B c → τ ν τ ) disfavours the scenarios with scalar operators. As can be seen across all the tables, our conclusions will not change much if we incorporate R J/ψ as input in our fit. This could be due to the large uncertainties present both in the predictions and measured values of R J/ψ .
Our goal in this analysis is to find out what type new operators can best explain the existing data. Our operator basis consists of all the linearly independent operators that are relevant for the b → cτ ν τ decays. Therefore, we did not extend our analysis to look for new physics models whose effects in b → cτ ν τ can be parametrised by one or more of the operators of our operator basis. Our priority is the analysis of the data, and the data can guide us to build models.
After ensuring that we are dealing only with the scenarios allowed by ∆AIC c , as well as constraints, we estimate the values of the parameters along with their uncertainties. Ideally, these would be obtained from the projections of the ∆χ 2 = 1 regions on the parameter line 3 . For simplicity, and to avoid asymmetric uncertainties, we consider a parabolic approximation around the chosen minimum and not only obtain the uncertainties of all parameters for each case, but also the correlation between them in the 2 parameter cases. These results are tabulated in tables 10, 11, 12 and 13. For some scenarios, instead of the results, the reader is asked to check the corresponding plots. In these scenarios, the best fit, instead of being an isolated point, is actually a contour in the parameterspace. Figure 1c is such an example. We note that we do not need large values of the C W s (< 1) to explain the observed discrepancies in general. Among the best possible scenarios, the data is more sensitive to the model with operator O V1 (with real C W ) or O V2 (with complex C W ) than the one with operator O T . From the best fit values we note that Re(C V2 ) < Re(C V1 ) << Re(C T ) < 1.
These results could be used by model builders to effectively put bounds on the parameter space of their lepton-flavour universality violating model, satisfying b → c ν transitions.

Prediction of observables and correlations amongst them
Using these NP results, we have predicted the values of the observables listed in table 1. Our predictions for all pertinent scenarios for the dataset without R J/Ψ are listed in tables 14 and 15. Predictions for the inclusive ratio R Xc are given in a separate table (table 16).
All of the predicted values for NP show deviations from their respective SM predictions. Moreover, neither all observables are equally deviated for a particular type of NP scenario, nor a single observable has similar deviations for different types of NP scenarios. Therefore, in trying to explain the deviations in R D ( * ) for a specific type of NP, we get information about the expected deviations in other associated observables. The obtained patterns then can be compared with the future measurements of these observables for a consistency check of the SM and to look for the types of NP. Any result, inconsistent with SM, but consistent with a future prediction of some observable, could be an indirect evidence in support for that specific scenario. In this regard, the correlations between the observables will play an important role. In figs. 3, 4, 5, and 6 we have shown the correlations between the various observables in different NP scenarios which are allowed by our model selection criteria. Following points illuminate our findings after scrutinising these plots: 1. Let us first note the very important correlations between R D and R D * in the scenarios which are allowed by our model selection criteria. In figs. 3a and 4a, we plot these correlations for the scenarios with one operator at a time, such as O T , O V1 , O V2 , and O S1 . We note that in all the scenarios, except the one with O V2 , the above two observables are positively correlated but the slopes are very different. By looking at the correlations, one can distinguish between the contributions from different NP operators. In the presence of either Re(C V1 ), Re(C V2 ), or Re(C S1 ), if one of the observables is consistent with the SM, then so must be the other. However, in the case of Re(C T ), there are regions in which R D * is consistent with the SM, whereas R D is largely deviated from its SM prediction. Therefore, if future data shows that R D * is within the SM ballpark but R D has a large value above its SM prediction, then any scenario with either Re(C V1 ), Re(C V2 ) or Re(C S1 ) has less chance to explain the data, but one with the operator O T will still be able to explain it. The situation is completely opposite in the case of Re(C V2 ), where the enhancement in R D * over its SM prediction is associated with a decrease in R D from its SM value. On the other hand, the contributions from O V2 with complex C V2 will show deviations in both the observables. 2. All the asymmetric and angular observables are insensitive to the operator O V1 , as its effect gets cancelled in    -Large deviations in both R D and R D * .
-P τ (D) will be consistent with its SM value.
-Measured values of P τ (D * ) and F L (D * ) are consistent with their respective SM predictions. -The measured value of A * F B will be above its SM prediction. 5. Let us accentuate a few other important points here.
In figs. 3b and 3c, we have shown the correlations between P τ (D ( * ) ) and R D * . In the presence of a tensor operator O T , the R D * is negatively and positively correlated with P τ (D) and P τ (D * ), respectively. However, when R D * is consistent with the SM, the τ polarization asymmetries will not be consistent with their respective SM predictions. In the presence of C T , the values of P τ (D) and P τ (D * ) will be below and above their respective SM predictions, respectively. Also, P τ (D * ) can be positive, whereas the SM predicted value is negative. In the same scenario, the correlations of R D * with A * F B and F L (D * ) are similar to those obtained for P τ (D * ) and P τ (D), respectively, for instances see figs. 3d and 3e. Also, here the forward-backward and the D * polarisation asymmetries are largely deviated from their respective SM predictions even when R D * is consistent with the SM. On the other hand, we do not see such behaviour in the presence of O S1 . In this case, the P τ (D ( * ) ), A * F B and F L (D * ) are consistent with their SM predictions depending on whether or not R D * is consistent with its SM prediction (see figs. 4b, 4c, 4d and 4e). 6. In case of the dataset with only R J/ψ dropped from the fit, the correlations of R D * with other observables like R J/ψ , R µ Λ , and A Λ F B are shown in figs. 3f, 3g,and 3h, respectively. Similar plots, which are obtained by dropping the BABAR data, are given in 4f, 4g, and 4h, respectively. In all the one-operator scenarios, the correlations are positive. Due to the large uncertainty in the SM prediction of R J/ψ , the predicted values of these observables are consistent with its SM prediction in all these cases. It is difficult to distinguish between the cases with either Re(C V1 ), Re(C V2 ) or Re(C S1 ). A very large deviation in R D * may allow us to see a small deviation in R J/ψ . Also, the contribution from Re(C T ) can be distinguished from other new operators. In a high precision experiment, contributions of various above mentioned operators are separable from each other by observing the correlation between R µ Λ and R D * . The contribution from Re(C T ) may allow a large deviation in R µ Λ , with a sizeable effect in R D * .
Similar patterns are observed in the correlations of A Λ F B . 7. Similar correlations in the allowed two-operator scenarios are shown in figs. 5 and 6. We note that it will be hard to distinguish the allowed two-operator scenarios from each other just from the correlations of R D * with R D , R J/ψ , and R µ Λ , as all the scenarios have similar correlations. However, the shape of the confidence regions the two-operator scenarios are different from those of the one-operator ones. , all other allowed twooperator scenarios are consistent with the SM even if there is a large deviation in R D * . By looking at these correlations, one will be able to distinguish a two-operator scenario from the one-operator ones. for two-operator scenarios previously showed in figures 3, 4, 5, and 6. We overlap these plots with the experimental results of RD * (eq. 1) and F D * L (eq. 22), upto 1σ (solid black) and 2σ (dashed black) ranges. The scenario with Re(CT ) is discarded in presence of the new result. [82]. According to that talk, the recently measured value is F D * L = 0.60 ± 0.08(stat.) ± 0.035(syst.) .
One immediate conclusion is that the operator O V1 alone can not explain such a large value in F L (D * ), since its effect is cancelled in the ratio. With the hope that this result may somewhat help us discriminate between our selected models, we recreate the F D * L vs. R D * correlation plots of figures 3, 4, 5, and 6 in figure 7. Keeping in mind the preliminary nature of this result, we show the corresponding experimental results up to 2σ. We note that the new tensor type operator with the Wilson coefficient C T cannot explain the observed result of D * polarisation asymmetry, though it is one of the best possible solutions for the explanations of the observed discrepancy in R(D * ). However, the other operators like O V2 alone, and the combinations of operators like O V1 , O S 1/2 are amongst the most probable scenarios that can accommodate the present data.

Summary
In this paper, we have predicted the SM values of the angular observables associated with the B → D ( * ) τ ν τ decays, following the results of an earlier up-to-date analysis on B → D ( * ) ν . Also, we have updated the SM prediction of R Xc using the results of [40] along with the proper correlations between the various non-perturbative parameters and masses. These predictions are based on two different schemes of the charm quark mass (M S and Kinetic). These include the NNLO perturbative corrections, and power-corrections up to order 1/m b 3 . We have separately mentioned results with power-corrections up to 1/m b 2 order as well. Our best results are R Xc = 0.214(4) for m c (3GeV ) = 0.987(13) GeV, and R Xc = 0.209 (4) when m kin c = 1.091 (20) GeV where in both the cases we have taken m kin b = 4.56 (21) GeV.
In the next part of our analysis, we have analysed the semitaunic b → cτ ν τ decays in a model independent framework with the ∆B = ∆C = 1 semileptonic operators. We have included the complete set of vector, scalar and tensor operators, while assuming the neutrinos to be left handed. Different possible combinations of all the effective operators have been considered, and following AIC c , the combinations, which are best suited for the available data, are considered for further studies. We have performed the analysis on several different prepared data sets. We note that for all of the data sets, the one-operator scenarios, with a real C W , can best explain the available data. However, in most of them, the scalar operators are not allowed by the constraint Br(B c → τ ν τ ) ≤ 30%. The most favoured scenarios are the ones with tensor (O T ) or (V − A) (O V1 ) type of operators. Also, the (V + A) type of interactions, with a complex C W , though less favoured, are allowed. In the absence of the BABAR data on R D ( * ) from our analysis, one-operator scenarios like (V ± A), S − P , and tensor operators with real C W are the most favoured ones. These one operator scenarios are easily distinguishable from each other by studying the correlations of R D * with R D and all the other asymmetric and angular observables. Also, the patterns of the future measurements of all such observables can easily discriminate the types of NP. Among all the possible combinations of (V ± A), tensor and (S − P ) operators, there are quite a few twooperator scenarios which pass all the selection criteria. In these cases, one cannot differentiate between the contributions from NP scenarios by looking at the correlations of R D * with R D , R J/ψ , and R µ Λ . However, the correlations of R D * with the various angular and asymmetric observables could be useful for such a discrimination. We have also predicted the numerical values of all the observables along with their errors, for the allowed scenarios.