Global fit to $b \to c \tau \nu$ transitions

We perform a general model-independent analysis of $b \to c \tau \bar{\nu}_\tau $ transitions, including measurements of $\mathcal{R}_D$, $\mathcal{R}_{D^*}$, their $q^2$ differential distributions, the recently measured longitudinal $D^*$ polarization $F_L^{D^*}$, and constraints from the $B_c \to \tau \bar{\nu}_\tau$ lifetime, each of which has significant impact on the fit. A global fit to a general set of Wilson coefficients of an effective low-energy Hamiltonian is presented, the solutions of which are interpreted in terms of hypothetical new-physics mediators. From the obtained results we predict selected $b \to c\tau\bar\nu_\tau$ observables, such as the baryonic transition $\Lambda_b \to \Lambda_c \tau \bar{\nu}_\tau$, the ratio $\mathcal{R}_{J/\psi}$, the forward-backward asymmetries ${\cal A}_\text{FB}^{D^{(*)}}$, the $\tau$ polarization asymmetries $\mathcal{P}_\tau^{D^{(*)}}$, and the longitudinal $D^*$ polarization fraction $F_L^{D^*}$. The latter shows presently a slight tension with any new-physics model, such that an improved measurement could have an important impact. We also discuss the potential change due the very recently announced preliminary $\mathcal{R}_{D^{(*)}}$ measurement by the Belle collaboration.


Introduction
The success of the Standard Model (SM) has reached its climax with the discovery of the Brout-Englert-Higgs boson [1][2][3], which seems to suggest the simplest scenario where the electroweak spontaneous symmetry breaking is linearly realized. In spite of its success as a low-energy effective field theory (EFT), there are both experimental signals and conceptual issues that cannot be accommodated in the SM framework and, therefore, motivate the search of New Physics (NP) beyond the SM. In this context, the series of anomalies in semi-leptonic B-meson decays, recently reported by several experiments, have caught a great attention in the scientific community. The unexpected deviations seem to appear in both b → c and b → s semi-leptonic decay transitions when different generations of leptons are involved, see Ref. [4] for a recent review.
The b → c transitions are of particular interest, because the necessary NP effect would be comparable with the tree-level contribution of the SM, which in turn would require NP to be either rather light or strongly coupled to the SM particles. Deviations from the SM predictions in those modes have been recently observed arXiv:1904.09311v1 [hep-ph] 19 Apr 2019 by the BaBar [5,6], Belle [7][8][9] and LHCb [10,11] collaborations in the ratios where B represents the branching ratio of the decay and denotes the light leptons, i.e., = e, µ. The combination of these measurements performed by the Heavy Flavour Averaging Group (HFLAV) [12] reads R avg D = 0.407 ± 0.039 ± 0.024 and R avg D * = 0.306 ± 0.013 ± 0.007 , with a correlation of −20%, which shows a tension of 4.4σ with our SM predictions (see also [12][13][14][15][16][17][18] to be discussed below. 1 Apart from the above observables, also the recent LHCb measurement [19] of the B c → J/Ψ ratio, deviates from the SM predictions R SM J/ψ ≈ 0.25-0.28 [20][21][22][23][24][25][26][27][28][29][30]. This points naively into the same direction, although the central value is in fact so large that it cannot be accommodated with NP contributions either. These deviations could be interpreted as hints at lepton flavour universality violation (LFUV), which cannot be accommodated in the SM and therefore suggest the existence of NP. The lack of evidence of similar discrepancies in K and π semi-leptonic and purely leptonic decays, or in electroweak precision observables, favours a scenario in which the potential NP contribution responsible for LFUV is only coupled to the third generation of leptons and quarks. The fact that in universality ratios large parts of the hadronic uncertainties cancel, renders underestimated theory uncertainties as an explanation extremely unlikely. However, recent measurements of R D * by LHCb [11] and Belle [9], which identify the final τ through its hadronic decays, result in values more compatible with the SM and yield a downward shift in the average that might suggest that the anomaly is smaller than indicated by the above numbers.
Our work aims at a better understanding of the nature of these anomalies. Instead of considering any specific NP model, we follow a bottom-up approach, in which the available experimental input is used to constrain any possible higher-scale effect and in this way infer information on NP without prejudice. We do exploit, however, the consequences of the apparent absence of NP close to the electroweak scale. Only afterwards we investigate which indications for more specific NP scenarios can be inferred. Numerous discussions can be found in the literature [16,, where the b → cτν τ transitions are studied from a model-independent point of view. However, most of these works restrict their analyses to either effects from a single NP operator or a single heavy particle mediating the interaction. We will adopt the most general possible scenario under a set of well-motivated assumptions instead.
In addition to the ratios defined in Eq. (1) we consider the normalized experimental distributions of Γ(B → D ( * ) τν τ ) measured by BaBar [6] and Belle [7]. Although this shape information was shown to provide quite stringent constraints in Ref. [6,37,41,44,60], it has been so far ignored in most phenomenological analyses. We also analyze the effect of including the recently announced value for F D * L by the Belle collaboration [61], F D * L = 0.60 ± 0.08 (stat) ± 0.04 (syst), which differs from its SM prediction by 1.6σ, and discuss its consequences in detail. Other related observables, such as P D * τ [9] and R J/ψ [19], are not included due to their large experimental uncertainties, 1 Note that this prediction does not rely on experimental inputs, but includes only part of the 1/m 2 q corrections in heavy quark effective theory. but are predicted from our fits. Very recently, the Belle collaboration has announced a new preliminary measurement of R D and R D * [62,63]: R Belle D = 0.307 ± 0.037 ± 0.016 and R Belle D * = 0.283 ± 0.018 ± 0.014 , with a correlation of −54%. This result is compatible with the SM at the 1.2σ level. Including this measurement in the global average yields which reduces the significance of the anomaly slightly; however, it still amounts to 4σ relative to the above SM prediction. 2 We present at the end of Sec. 3 an updated analysis, including in the fit these preliminary data, and discuss their implications.
Our paper is organized as follows: In Sec. 2, the theoretical framework used in this work is presented, and the physical observables and experimental inputs are defined. In Sec. 3, we discuss our global χ 2 fit and detail the resulting values of the fitted parameters. The interpretation of these results and their relation to NP are given in Sec. 4, where we complete our discussion with several additional fits, relaxing some of the assumptions. A set of predictions for relevant observables, for which measurements will be published or improved soon, is presented in Sec. 5. Finally, we draw conclusions in Sec. 6. Some technical details are relegated to the appendices.
The above fermionic operators are given by 3 and are weighted by the corresponding Wilson coefficients C i , which are, in general, lepton and flavour dependent, and parametrize any possible deviation from the SM, i.e., C SM i ≡ 0. This effective Hamiltonian forms the basis of our analysis, restricted only by a minimal set of well-motivated assumptions: • Possible NP contributions are assumed to be present only in the third generation of leptons. This is motivated by the absence of experimental evidence of deviations from the SM in tree-level transitions involving light leptons; specifically, precision measurements like the ratio B(τ → µν τνµ )/B(τ → eν τνe ) = 0.9762 ± 0.0028 [64] and the analysis of b → c(e, µ)ν (e,µ) transitions in Ref. [52] constrain potential effects to be negligible in the present context.
• The coefficient C V R is assumed to be lepton-flavour universal in our main fit. This statement can be derived [65][66][67] in the context of the Standard Model Effective Field Theory (SMEFT) [68,69], which is the appropriate effective theory in the presence of a sizable energy gap above the electroweak scale if the electroweak symmetry breaking is linearly realized. The experimental facts that no new states beyond the SM have been found so far up to an energy scale of approximately 1 TeV and that measurements of the Higgs couplings are all consistent with the SM expectations support this scenario. In this case, C V R is strongly constrained from b → c(e, µ)ν (e,µ) data [52], and we set it to zero for convenience. If the assumption of linearity is relaxed, a non-universal C V R coefficient can be generated [67]; we will consider this case separately.
• The CP-conserving limit is taken, so all Wilson coefficients C i are assumed to be real. This is mostly done for convenience; however, none of the measurements related to the B anomalies refers to a CPviolating observable. Possible CP-violating contributions have been analyzed before in, e.g., Ref. [44,58,59,70,71]. Note that in the presence of such couplings other observables can become relevant, like electric dipole moments, see, e.g., [72,73]. This assumption will be briefly commented in Section 3.

Form Factors
The relevance of hadronic uncertainties in the determination of |V cb | has opened an intense debate about the most adequate way to parametrise the relevant hadronic form factors [15][16][17][18][74][75][76][77]. It has been suggested that the accuracy of the usually adopted Caprini-Lellouch-Neubert (CLN) parametrisation [78] has been probably overestimated and the current experimental precision requires to use more generic functional forms such as the one advocated by Boyd, Grinstein and Lebed (BGL) [79][80][81]. However, we note that the observables considered here are mostly ratios, reducing the overall form-factor sensitivity. We consider a heavy quark effective theory (HQET) [82,83] parametrization, including corrections of order α s , Λ QCD /m b,c and partly Λ 2 QCD /m 2 c , mostly following [16,52]. In the heavy-quark limit all form factors either vanish or reduce to a common functional form, the Isgur-Wise function ξ(q 2 ) [84]. Thus, it is convenient to factor out ξ(q 2 ) by defining [16]ĥ (q 2 ) = h(q 2 )/ξ(q 2 ) .
The leading Isgur-Wise function can be more conveniently expressed in terms of the kinematical parameters The variable ω(q 2 ) is the inner product of the B and D ( * ) velocities, so that ω = 1 corresponds to the zero-recoil point, q 2 max = (m B − m D ( * ) ) 2 , where ξ(q 2 max ) = 1. The conformal mapping z(q 2 ) encodes in a very efficient way the analyticity properties of the form factors, transforming the cut q 2 plane into the circle |z| < 1 [85], so that a perturbative expansion in powers of z(q 2 ) has an optimized convergence. Up to O(z 4 ) corrections, ξ(q 2 ) can be written as 4 (12) and it is characterized through the parameters ρ 2 , c and d.
The functionsĥ(q 2 ) introduce corrections of order Λ QCD /m b,c and Λ 2 QCD /m 2 c via the subleading Isgur-Wise functions χ 2,3 (ω), η(ω) at order 1/m c,b and l 1,2 (ω) at order 1/m 2 c , parametrized by the parameters {χ 2 (1), χ 2 (1), χ 3 (1), η(1), η (1)} and {l 1 (1), l 2 (1)}, respectively. They also include the corrections of order α s . The detailed parametrization of the different form factors can be found in Ref. [16,52]. The main difference to the latter article is the introduction of the z 3 term in the leading Isgur-Wise function, that renders the fit compatible with the extrapolation of the recent lattice data [74,87] to large recoil. We updated the corresponding fit to the inputs from lattice quantum chromodynamics (QCD) [74,[87][88][89], light-cone sum rules [90] and QCD sum rules [91][92][93] (see [52] for details); note that this fit does not make use of experimental data, thereby rendering the form factors independent of the NP scenario considered. The results obtained for the 10 form-factor parameters are given in Table 1, while the corresponding correlation matrix can be found in Table 10 of Appendix A.

Parameter
Value 0.14 ± 0.23 l 2 (1) 2.00 ± 0.30 Table 1: Inputs used to determine the form factors in the HQET parametrization as in [16]. The first three parameters determine the leading Isgur-Wise function, while the last seven enter in the 1/m c,b and 1/m 2 c corrections. The correlations between these parameters can be found in Table 10 of Appendix A.

Observables and experimental input
We collect the formulae for the main observables entering our analysis. Starting with B → D ( * ) τν τ decays, we obtain from the effective Hamiltonian of Eq. (8) their differential decay rates as a function of the general set of Wilson coefficients [35,94]: The helicity amplitudes, which encode the information from the hadronic form factors, can be found in the Appendix B. The values of the quark and meson masses and other experimental inputs used in our analysis are listed in Table 9 of Appendix A.
Besides the semi-leptonic processes included in the fit, the pure leptonic decay B c → τν τ is crucial in determining the direction of potential NP effects, since it strongly constrains the axial (C V R − (1 + C V L )) and, especially, the pseudo-scalar (C S R − C S L ) contributions [42,95]: From these expressions, four classes of observables are obtained that are determined in experimental analyses: • The ratios R D ( * ) Experimental measurements of the ratios R D and R D * have been published by BaBar [5,6], LHCb [10,11], and Belle [7][8][9] (see also [62]) using different techniques. These results have been averaged by the HFLAV collaboration, giving the values listed in Eq. (2) [12]. The results for each experiment and their average are also shown in Fig. 1, with and without the result from Ref. [62,63].
As mentioned above, these ratios are advantageous both theoretically and experimentally, as they allow for the cancellation of uncertainties, specifically the CKM factors and leading form factor uncertainties on the theoretical side.
• Differential distributions of the decay rates Γ(B → D (*) τν τ ) Belle and BaBar have also provided data on the measured q 2 distributions for B → D ( * ) τν τ [5,7]. We show the reported binned values in Appendix A, Table 8. Since the global normalizations of these distributions are effectively already included via the values for R D ( * ) in these analyses, they are not independent degrees of freedom. This can be taken into account either by introducing a free normalization factor for the distributions as in Ref. [44] or by normalizing the differential binned distributions in the following way: which keeps the information about the shape of the distribution, independently of the global normalization. The treatment of systematic uncertainties and correlations follows Ref. [44].
• The leptonic decay rate B c → τν τ While this decay is not expected to be measured in the forseeable future, it can still be used as a constraint in the following way: A 30-40% upper bound can be derived from the B c lifetime [42,44,96]. A more stringent 10% bound has been recently obtained from LEP data at the Z peak [97], and it may become even stronger by performing the analysis with the full L3 data [98]. However, this bound assumes the probability of a b quark hadronising into a B c meson to be the same at LEP (e + e − ), the Tevatron (pp) and LHCb (pp), which exhibit very different transverse momenta. This is known to be a bad approximation in the case of b-baryons, see Ref. [12]. The bound also makes use of the SM theoretical prediction for B(B c → J/Ψ ν ). See also Ref. [58] for a more detailed discussion.
In our fits, we will compare the two options of imposing the upper bounds B(B c → τν τ ) < 10% (30%). The bounds are used in a way that only points in the parameter space that fulfill this constraint will be considered. 5 • The longitudinal polarization fraction F D * L A measurement of the D * longitudinal polarization fraction, defined as has been recently announced by the Belle collaboration [61]. The explicit expression for Γ λ D * =0 (B → D ( * ) τν τ ) is given in Appendix C. Being normalized to the total rate, this observable also enjoys the advantages of the other ratios. To study the implications of this measurement, we perform one fit with it and one without it.

Fit and results
In order to extract the information on the NP parameters C i , we perform a standard χ 2 fit. The χ 2 function can be splitted in two parts, where χ 2 exp contains the experimental information discussed in the last subsection (again a sum of the three main contributions) and χ 2 FF the information on the form factors discussed in Sec. 2.2 in the form of pseudoobservables with the "experimental" information presented in Table 1. Each individual χ 2 is defined as: with y i denoting the input parameters of the fit, i.e., , l 2 (1)}, ρ ij the correlation between the observables i and j, and σ i the uncertainty of the observable i. In the above equation, f th represents the theoretical expression for a certain observable and f exp its experimental value. The contribution from the limit on the branching fraction of B c → τν τ is implemented as a Heavyside Theta function, its contribution being zero for parameter combinations where the limit is obeyed and infinity for those where it is not. The uncertainty of a parameter y i is determined as the shift ∆y i in that parameter, where the minimization of χ 2 | y i =y min i +∆y i varying all remaining parameters in the vicinity of the minimum leads to an increase of ∆χ 2 = 1.

Standard Model
We start by discussing the situation in the SM, corresponding to C i ≡ 0. The global fit to the data discussed above does actually appear to be reasonable: we obtain χ 2 min = 65.5 for, naively, 57 degrees of freedom (d.o.f.), corresponding to a naive confidence level (CL) of ∼ 20%. However, these numbers are misleading for the following reason: The systematic uncertainties added to the dΓ/dq 2 distributions have been chosen to be maximally conservative. Therefore, it can be expected that the corresponding χ 2 contribution is reduced; this is indeed seen since the contribution from these distributions is χ 2 min,dΓ ∼ 43 for, again naively, 54 d.o.f.. Considering instead the contribution from R D ( * ) , we do of course reproduce the well-known puzzle, i.e., we obtain χ 2 = 22.6 for 2 d.o.f., corresponding to a 4.4σ tension. Note also that the limit from the B c lifetime is irrelevant in the SM fit.
These observations imply that also NP scenarios should not be judged simply by χ 2 vs. d.o.f., but by the improvement they yield when compared to the SM.

New Physics
Since the Wilson coefficients enter each observable bilinearly (the coefficient of the left-handed vector operator being (1 + C V L )), there is a degeneracy between a set of Wilson coefficients and a mirror minimum with The two sets of Wilson coefficients give identical predictions for all observables and consequently have the same χ 2 value. 6 In the following, we will always discuss the closest minimum to the SM scenario, i.e., with smaller |C V L |, and will omit the sign-flipped solution; this corresponds to considering only values The global fit to the data described in Section 2.3 without including the longitudinal polarization yields a unique global minimum (for C V L > −1) with χ 2 Min 1 = 34.1 for 53 d.o.f.; in addition, we find two local minima, with χ 2 Min 2 = 37.5 and χ 2 Min 3 = 58.6, the latter of which is, however, highly disfavoured by the differential distributions. We summarize the results for the NP parameters in Table 2. Including the recently announced longitudinal polarization in the global fit, we find that the overall structure for the lower two minima remains the same; however, this observable reduces slightly the available parameter space for the NP parameters. The central values of the scalar NP parameters are smaller for the global minimum, while the 1σ-ranges remain almost constant. The most striking effect is that the already less favoured local minimum disappears. The results for the NP parameters in this context can be found in Table 3. In both cases the form factor parameters reproduce their input distributions up to very small shifts. For illustration we show graphically in Fig. 2 the NP parameters for the different minima obtained in the two scenarios. There are important correlations between the NP parameters obtained from the fit. We illustrate them in the two-dimensional plots in Fig. 3 for the different scenarios. The contours shown there are relative to the global minimum.
We note that the distributions for, especially, the scalar parameters are highly non-gaussian. Reasons are the way the upper limit on B(B c → τν τ ) is included and the fact that the first two minima overlap to some extent. The former is also the reason for the strong asymmetry in the uncertainties for C S L,R . Since only their sum and difference enter B → D and B → D * decays, respectively, these parameters are furthermore highly correlated. The local minima are not very deep, resulting in complications in the determination of   Table 2: NP parameters for the minima obtained from the χ 2 minimization and 1σ uncertainties. There are, in addition, three corresponding sign-flipped minima, as indicated in Eq. (20). In the first three columns, the constraint B(B c → τν τ ) ≤ 10% has been applied, whereas in the last three columns, this requirement has been relaxed to B(B c → τν τ ) ≤ 30%.
the uncertainties for the Wilson coefficients at these points.

Min 1b Min 2b
Min 1b Min 2b  Table 3: NP parameters for the minima obtained from the χ 2 minimization including F D * L and their 1σ uncertainties. There are, in addition, the corresponding sign-flipped minima, as indicated in Eq. (20).
The fit results for the R D and R D * ratios at the different minima are presented in Fig. 4. As expected, the predictions obtained from the fit are compatible at the 1σ level with the experimental data, in the case of Min 1 and Min 1b essentially reproducing them. From the fit results without including F D * L , the following information can be extracted: • The reduction of the global χ 2 by 31.4 (31.7) for 4 NP parameters implies a strong preference of NP compared to the SM, taking the present data set at face value and B(B c → τν τ ) ≤ 10% (30%).
• There is no absolute preference of a single Wilson coefficient in the sense that for the global minimum each individual Wilson coefficient is compatible with zero within at most 1.1σ.
• On the other hand, considering scenarios with only a single Wilson coefficient present, there is a clear preference for C V L : removing the other three Wilson coefficients increases χ 2 only by 1.4, corresponding to 0.14σ. Hence, Min 1 is well compatible with a global modification of the SM, that is, C V L being the only non-zero coefficient. • The other two minima are numerically further away from the SM; instead of a single dominant contribution, there are several sizable Wilson coefficients whose contributions partly cancel each other in some observables. These minima also imply different values for the fitted observables: Min 2 corresponds to a slightly worse fit for both, R D ( * ) and their q 2 distributions, while Min 3 fits R D ( * ) perfectly, but is essentially already excluded by the (rather coarse) measurements of the distributions available.
• All minima saturate the constraint B(B c → τν τ ) ≤ 10% (30%). Relaxing the upper bound allows for a larger splitting between the two scalar Wilson coefficients, and the contribution of the scalar operators gets enlarged. This constraint is consequently the main argument at low energies disfavouring a solution with only scalar coefficients. Any such solution would require a lower value for R D * by about 2σ.
• Having solutions with relevant contributions from all Wilson coefficients illustrates the importance of taking into account scalar and tensor operators in the fit.
• The fit results for the form factor parameters reproduce their input values displayed in Table 1 up to tiny shifts. This implies that the uncertainties of the experimental data with tauonic final states are large compared to the hadronic uncertainties. Differently stated, while the ranges obtained for the NP parameters are obtained in fits varying all form factor parameters simultaneously with the NP ones, they are essentially determined by the experimental uncertainties at the moment.
• Generalizing the fit to complex Wilson coefficients does not improve the minimal χ 2 value, but opens up a continuum of solutions. Hence complex Wilson coefficients can explain the anomalies as well as real ones, but they do not offer any clear advantages regarding the fit quality, so they have not been considered here for simplicity. It should be mentioned, however, that in specific models the option of complex Wilson coefficients can open up qualitatively new solutions, as for example the model proposed in Ref. [70], where only the coefficients C S L ,T (C S L ∼ C T ) are present, requiring a non-vanishing imaginary part in order to accommodate the experimental data. This fact implies correlations with new observables like electric dipole moments, which can then be used to differentiate this model from solutions allowing for real coefficients [73].
• As discussed above, for each minimum given in Table 2 there is a degenerate solution, see Eq. (20).
Including the recent measurement of the longitudinal polarization F D * L in the global fit, the above statements hold up to the following differences: • Still there is no clear preference for a single Wilson coefficient. The central values for the scalar coefficients are smaller for the global minimum, such that the bound from the B c lifetime is not saturated even in the 10% case. As a consequence, the minimum does not change when allowing for larger values of B(B c → τν τ ), only the allowed parameter ranges increase.
• The second local minimum (previously referred to as Min 3) disappears.
It is not straightforward to compare our fit with the results from other analyses in the literature, because we are including the information from the q 2 distributions that has been ignored in previous fits with the exception of Ref. [37,41,44,60]. Besides that, some works include additional observables such as R J/ψ or slightly different bounds on B(B c → τν τ ). Nevertheless, comparing the findings of previous fits with our results is quite enlightening since it illustrates the relevance of the additional observables we are considering. Generic fits to the R D ( * ) world averages in Eq. (2), with the effective Hamiltonian of Eq. (8) [16,31-56, 58,59], have shown the existence of many possible solutions, some of them involving only one or two Wilson coefficients. Including the B(B c → τν τ ) upper bound reduces the number of allowed possibilities, but several different scenarios remain still consistent with the data. Dropping the binned q 2 distributions from our fit, we can easily reproduce all those solutions. However, most of them lead to differential distributions in clear conflict with the BaBar and Belle measurements. While a sizable new-physics contribution to some Wilson coefficient can easily generate the needed enhancement of the B → D ( * ) τν τ rates, it tends to distort the shape of the differential distributions in a way than can no-longer accommodate the data, similarly to what happens for Min 3. Once the full experimental information on R D ( * ) (rates and binned distributions) is taken into account, the χ 2 minimization only gives the three solutions shown in Table 2, and when including F D * L in the fit, the number of solutions is further reduced to two.
Finally, a few comments on the very recent measurement of R D and R D * released in Moriond by Belle [62,63] are in order. It should be kept in mind, however, that these results are still preliminary. Including the new average in the fit (see Fig. 1), we find again qualitatively similar solutions as before, as can be seen by comparing the numerical results in Tables 3 and 4. We show for simplicity only the solutions with B(B c → τν τ ) < 10%; increasing this limit results again essentially in larger ranges for especially the scalar Wilson coefficients, although the new global minimum now does saturate this limit, so also the central values do change. Again all individual coefficients are roughly compatible with zero at 1σ. C V L alone also still provides an excellent fit to all the data, now with a smaller central value of ∼ 0.08. Interestingly, the fit with only C T is improved by the new results, which, however, does not correspond to a simple singlemediator scenario, as discussed below. However, related to that observation, also the fit in the scenario of Ref. [70] improves by ∆χ 2 = −1.8 (for B(B c → τν τ ) < 30%).

Interpretation of results
In Sec. 3 we have described the global fit to the available data on b → cτν τ transitions in terms of the Wilson coefficients of an EFT framework defined at the b-quark mass scale. The EFT in this range is conventionally  Table 4: Minima and 1σ uncertainties obtained from the global χ 2 minimization, including the new preliminary result measured by Belle on the R D ( * ) ratios and the F D * L polarization, using B(B c → τν τ ) < 10%. There are, in addition, the corresponding sign-flipped minima, as indicated in Eq. (20). called Weak Effective Theory (WET) and is composed of the five lightest quarks and the three generations of leptons, and ruled by the SU (3) C ⊗ U (1) Q gauge symmetry. This is a valid approach assuming -as strongly suggested by all available collider data -that no new degree of freedom exists coupling to this channel with a mass around or lower than the b quark. However, ultimately the goal is to gain insight into the high-energy structure of the theory. To that aim, renormalization-group techniques are used to relate the coefficients extracted in our analysis to those relevant at the scale of the potential new high-energy degree(s) of freedom. This process involves several scales and thresholds, see Fig. 5.
The relation to the coefficients at the electroweak scale is determined by QCD and are known [99][100][101][102]. Above the electroweak (EW) scale, the Lagrangian has not undergone spontaneous symmetry breaking and, therefore, the fermionic fields should be expressed in terms of weak eigenstates rather than mass eigenstates. Moreover, the top quark, the electroweak gauge bosons and the Higgs boson have to be considered as new degrees of freedom in the theory. The relevant framework at this scale is the full SM, with the addition of the effects of NP. For relatively low NP scales 1 TeV, the relevant new degrees of freedom can be included explicitly. However, the suggested absence of new degrees of freedom below ∼ 1 TeV allows us to parametrize any NP contribution in the framework of another effective theory. This can be the so-called SMEFT under the conditions specified in Section 2, or a more general framework with a non-linear representation for the Higgs, see, e.g., Ref. [103,104].
In SMEFT, the effective lagrangian can be expanded in inverse powers of the NP scale, Λ NP , i.e., built from a series of higher-dimensional operators in terms of the SM fields and invariant under the SM gauge group SU (3) C ⊗ SU (2) L ⊗ U (1) Y [68]. A convenient complete and non-redundant basis of dimension-six operators is the Warsaw basis [69]. In order to relate both EFTs, the matching between the WET theory and the SMEFT has to be performed at the EW scale [65,66,101,102,105,106]. The matching onto the basis in the non-linear case [107,108] is given in Ref. [67].
As an illustration of the effect of the running, we show the relation between the WET Wilson coefficients at µ b ≈ 5 GeV and the SMEFT Wilson coefficients at an hypothetical NP scale of Λ = 1 TeV, calculated in Ref. [55,112], which can be trivially inverted: For a discussion of the notation used for the SMEFT Wilson coefficients in the Warsaw basis see Appendix E.
With the coefficients at the potential NP scale at hand, one can try to go beyond the EFT framework and get an idea about which class of NP might be responsible for the observed pattern: at the scale Λ, the coefficients C i should result from integrating out the new heavy degrees of freedom. In Table 5, the quantum numbers of all possible candidates able to participate in the b → c transitions are listed and their nature is identified (see also [37]). We note that, in some cases, a given NP mediator may contribute to more than one Wilson coefficient, thus resulting in correlations among them. In Appendix D, we list the effective Lagrangians obtained after integrating out each of the possible heavy degrees of freedom. We show in the last two columns of Table 5 the set of Wilson coefficients to which the new degrees of freedom contribute, both in the SMEFT and in the WET. The RGE running changes the relative size of these coefficients, as seen above, and causes mixing among the operators O S L and O T . When considering such specific classes of models, generally other constraints apply. Specifically, searches for the corresponding mediators can exclude a large part of the parameter space, or even the whole scenario (like the W ) [114][115][116]. In the following we will not discuss these constraints, but simply give examples for how the required coefficients could be generated, irrespective of their actual viability.
We are now in a position to interpret the different solutions obtained in the fit shown in Table 2 and  Table 3. Let us focus first on the scenarios where F D * L is not included. The minimum with highest χ 2 , Min 3, presents relevant contributions from the operators O S L and O T . The origin of these Wilson coefficients could be explained, for instance, with the presence of the scalar leptoquarks R 2 ∼ (3, 2, 7/6) or S 1 ∼ (3, 1, 1/3), whose contributions to the Lagrangian at the NP scale are given in Appendix D. An additional mediator would be necessary to generate the sizeable contribution to C V L , however, in the former case. Min 2, which exhibits non-zero values for all Wilson coefficients, could be explained by combinations of several candidates, for instance S 1 and H 2 . Also for Min 1 there are different possibilities, since the fit Spin Q.N. Nature Allowed couplings SMEFT WET  does not single out a specific coefficient. However, the simplest option remains the scenario where the only relevant contribution is proportional to the SM one, i.e., all Wilson coefficients but C V L are compatible with zero at 1.1σ. This possibility could be generated, for instance, by the effect of a W boson, with M W /(g ν g † du ) 1/2 ∼ 2 TeV. For a sequential W with SM couplings, one would need M W ∼ 0.2 TeV, which is already ruled out by direct searches [117]. More exotically, but more realistically given the aforementioned high-energy constraints, one could explain the modification on the O V L operator by introducing leptoquarks (LQs), such as the vector U 3 ∼ (3, 3, 2/3) or the scalar S 1 ∼ (3, 3, 1/3) LQs. However, extra symmetries in the UV regime would have to be assumed in order to guarantee that other flavour transitions compatible with the SM are respected.
In Fig. 6 we show the dependence of selected observables on individual Wilson coefficients. The left-top panel in Fig. 6 shows that it is straightforward to achieve consistency with the experimental measurements for R D ( * ) by shifting only the Wilson coefficient C V L , i.e., modifying the SM coefficient. The polarization observables show a good potential to differentiate between different contributions. Particularly interesting is the longitudinal polarization fraction in B → D * τν τ , shown in the bottom-right panel, for which the Belle collaboration recently announced a first measurement [61]. As this sub-figure shows, it is difficult to accommodate it at 1σ for any of the individual Wilson coefficients. The only contributions allowing for a significantly larger value of this observable than in the SM are those from scalar operators; however, values accommodating F D * L are in conflict with the bound from B(B c → τν τ ) < 10% (dashed lines), and extending this bound to 30% still does not allow to accommodate its central value. This figure therefore indicates why none of the fit scenarios yields values for F D * L in the 1σ range; we take this as a motivation to investigate the consistency of the different measurement in more detail.
In order to do so, we use the fact that only three combinations of the four Wilson coefficients enter B → D * τν τ observables as well as the leptonic B c decay: C V L , C T and the pseudo-scalar coefficient C P ≡ C S R − C S L . Every observable therefore results in a non-trivial constraint in the C P − C V L plane if C T is fixed to some value. We show the preferred parameter ranges obtained for the individual observables  in Fig. 7, for a representative set of C T values. The combination of R D * and the bound on B(B c → τν τ ) determines a narrow strip in this parameter plane, dominated by the former for the bound on C V L and the latter for the bound on C P . The overlap of the other observables varies with the value for C T ; however, there is no value of C T for which all 1σ bands overlap. In fact, the 1σ range for F D * L cannot be reached by any NP parameter combination in this setup, when only imposing the B(B c → τν τ ) constraint of 10% or even 30% and at the same time requiring a positive shift in R D * . Agreement can presently be achieved at the 2σ level; nevertheless, a confirmation of the present central values with higher precision could indicate the inconsistency between the data and any NP with flavour-universal C V R .
This potential incompatibility would suggest one of several possibilities: 1) One of our theoretical assumptions is incorrect and the SMEFT cannot be applied at the electroweak scale. This could happen if one or several of the following cases apply: (a) There is an insufficient gap between the electroweak and the NP scale, i.e., there are new degrees of freedom close enough to the EW scale to invalidate an EFT approach. (b) The electroweak symmetry breaking is non-linear, changing also the character of the observed Higgs-like particle. In that case C V R could contribute to the fitted observables, because it would no-longer be necessarily flavour universal. (c) There are additional light degrees of freedom like right-handed neutrinos [118][119][120], yielding additional operators.
Note that we also assumed the semi-leptonic decays with light leptons to be free from NP. However, the corresponding constraints are so strong that even relaxing this assumption would not significantly change our analysis [52].
2) An unidentified or underestimated systematic uncertainty in one or several of the experimental measurements.
In any case, the upcoming experimental studies of not only the LHCb collaboration, but also the Belle II experiment which started to take data will hopefully resolve this question soon. For completeness of our discussion, we have consequently performed the fit relaxing the condition of flavour universality on C V R . As a consequence of adding C V R as an extra d.o.f. to fit, the number of solutions is enlarged. As shown in Fig. 8, one finds now four different solutions (plus their sign-flipped counterparts), given numerically in Table 6 Table 6: Minima with their 1σ uncertainties obtained from the global χ 2 minimization, including F D * L and B(B c → τν τ ) < 10% in the fit while allowing for C V R = 0. There are, in addition, the corresponding sign-flipped minima, as indicated in Eq. (20).
The doubling of minima can be understood qualitatively in the following way: B → D is dominated by the combination of Wilson coefficients corresponding to the vector coupling C V = 1 + C V L + C V R , while B → D * is dominated by the axial-vector coupling C A = C V R −(1+C V L ). Their rates are correspondingly roughly given by |C V,A | 2 . For C V R ≡ 0 we have C V = −C A , and the only remaining discrete symmetry is that discussed in Section 3.2, the second solution being eliminated by our choice C V L > −1. With a finite coefficient C V R , these two solutions become four ({C A = ±|C A |, C V = ±|C V |}), since now |C A | = |C V |; two of those are again eliminated by our choice for C V L , leaving two solutions per minimum with C V R ≡ 0. This degeneracy is broken by interference terms, notably Re(C A C * V ) in B → D * , but also the interference with scalar and tensor operators. Nevertheless, this approximate degeneracy explains the doubling of solutions for finite C V R .
As can be seen from the comparison of Table 6 with Table 3, the previous global minimum, Min 1b, remains a solution of this more general fit, now called Min 6. Min 7 is again relatively close to Min 6, however with a significant contribution from C V R and hence qualitatively different from Min 2 in the previous fits. The new global minimum Min 4 and the close-lying Min 5 improve the agreement of the fit with the data significantly. However, in these scenarios the SM coefficient is almost completely cancelled and its effect replaced by several NP contributions. These are hence fine-tuned scenarios, and should be taken with a grain of salt. We have also analyzed the individual observables in B → D * and the bound on B(B c → τν τ ) for this case. This is illustrated in Fig. 9, for different benchmark values of C V L and C T , in the plane C V R − C P . The figure shows again the allowed regions at 1σ for the different observables. In accordance with the above reduction for χ 2 min , we observe that in this case it is possible to have an overlap of all the bands. However, it is still not possible to reach the central value for the longitudinal polarization fraction, and as mentioned above, this scenario corresponds to a highly fine-tuned combination of parameters.

Predictions
We use our global fits from Sec. 3 to predict selected observables that are either not measured yet, but expected to be measured soon, or presently measured with uncertainties that are larger than those from the fits. These additional measurements serve two purposes: firstly, they provide additional information that is theoretically related, but experimentally independent (to varying extent) from existing measurements, thereby helping to establish NP and excluding underestimated systematic uncertainties as the source for the anomaly. Secondly, they can provide experimental information on combinations of Wilson coefficients that are not or only weakly constrained so far, thereby allowing to distinguish different NP scenarios.
We will first present the predictions for observables of the key modes B → D ( * ) τν τ , before focusing on other semi-leptonic decays, specifically Λ b → Λ c τν τ and B c → J/ψ τν τ .

Predictions for B → D ( * ) τν τ observables
We start by analyzing the q 2 distributions of several angular observables. While these distributions can be very effective in distinguishing different NP scenarios, they are difficult to measure, due to the missing information on the neutrinos. The angular dependence of the differential decay width B → D ( * ) ν can be parametrized by three independent angular coefficients, which are in principle experimentally accessible. Here, θ is the angle between the D ( * ) and chargedlepton three-momenta in the -ν center-of-mass frame. An angular observable commonly defined in the literature is the forward-backward asymmetry, which is determined by the b ( * ) (q 2 ) coefficient according to the following expression: This observable yields complementary information, since it does not contribute for quantities integrated over the full range of cos θ . One can also decompose the differential branching ratio according to the two possible polarizations of the charged (τ ) lepton, giving rise to another observable named τ polarization asymmetry: where λ τ is the helicity of the τ lepton, and dΓ D ( * ) λτ /dq 2 is the differential decay width of B → D ( * ) τν τ for a given helicity λ τ .
Analogously, one can extract from the angular distribution in the secondary D * → Dπ decay the fraction of longitudinally polarised D * mesons by constructing the following observable: In Fig. 10, we show the q 2 dependence of the B → D ( * ) τν observables defined above, for the two solutions obtained in the global fit including F D * L , Min 1b and Min 2b, together with their SM prediction. Using these observables, Min 2b could rather clearly be differentiated from both the SM and Min 1b. The same is not true for Min 1b and the SM, for the simple reason that this minimum is compatible with only shifting the SM coefficient at 1σ. In that case the SM predictions are unchanged, which means that the width of the red bands is due to the possible presence of additional NP operators. Precise measurements of these distributions could hence show the existence of operators other than O V L .
Given the aforementioned difficulty with measuring q 2 distributions, typically the integrated observables are measured first, defined as where O(q 2 ) refers to the numerator in the ratios, i.e., numerator and denominator have to be integrated separately. The Belle collaboration has in fact released results for two integrated quantities, the τ polarisation asymmetry P D * τ = −0.38 ± 0.51 (stat) +0.21 −0.16 (syst) [121], and the recently announced longitudinal polarisation of the D * meson, F D * L = 0.60 ± 0.08 (stat) ± 0.04 (syst) [61,122]. In Fig. 11, we show the predictions for the integrated observables of B → D ( * ) τν τ , together with their experimental values where available. Clearly already the integrated observables provide a possibility to distinguish the different NP scenarios. The fitted values for F D * L are closer to the experimental results for the fits including this observable, which is to be expected. However, they fail to reproduce the measurement within 1σ, as discussed above, which renders a more precise measurement of this quantity an exciting prospect.

Predictions for R Λc
Another observable that could shed light on the R ( * ) D puzzle is the Λ b → Λ c τν τ decay, in particular the universality ratio This decay mode has not been observed yet, but LHCb has the potential to perform this measurement in the near future.
On the theoretical side, the differential decay rate Λ b → Λ c ν has been calculated in terms of the helicity amplitudes [123,124]: where The superindices V A indicate vector and axial-vector contributions (C V R ± C V L ), SP scalar and pseudoscalar (C S R ± C S L ), and T tensor contributions (C T ). Being a baryonic decay, this mode is sensitive to different combinations of Wilson coefficients than B → D ( * ) τν τ . We use the parametrization of the QCD form factors from Ref. [123,124], which take the simple form: The numerical values of the corresponding form-factor parameters, extracted from lattice data [123,124], are displayed in Table 7. Other relevant experimental inputs are summarized in Table 9. Fig. 12 shows the predicted ratio R Λc and its uncertainty for the three minima of Table 2 (Min 1, Min 2 and Min 3), the two minima including F D * L of Table 3 (Min 1b and Min 2b), with the upper limit B(B c → τ ν) ≤ 10% and the SM prediction. The errors considered here just take into account the variation of the Wilson coefficients and the parametric error for the lattice input. Other systematic errors are not shown.   Figure 12: Predictions for R Λc (left) and R J/ψ (right) for the minima of Table 2 and Table 3, with an upper bound B(B c → τ ν) ≤ 10%. The SM prediction is shown as a blue band. The experimental value of R J/ψ is given by the gray band.

Predictions for R J/ψ
The ratio has been recently measured by LHCb with the run-1 dataset (3fb −1 ) [19]. We have not included this observable in our fit because the hadronic uncertainties are not at the same level as for the observables related to B → D ( * ) transitions and the experimental error is large. Instead, the predictions for this observable are computed and compared with the current data. The experimental uncertainties are expected to be significantly reduced with the larger statistics already accumulated at LHCb. The differential decay rate for this transition can be expressed in a similar way than theB → D * distribution in Eq. (14) [125]: where λ J/ψ (q 2 ) = (m Bc − m J/ψ ) 2 − q 2 (m Bc + m J/ψ ) 2 − q 2 is the usual Källén function and H i are the hadronic helicity amplitudes, similar to the ones used for the decay rates of Sec. 2, which can be found in Appendix B.
The predicted values of R J/ψ for the minima of Tables 2 and 3 as well as for the SM, are given in the right panel of Fig. 12. Again the errors considered here just take into account the variation of the Wilson coefficients and the parametric error for the lattice input. For this observable, there are additional theoretical uncertainties associated with the parametrization of the form factors, which are difficult to quantify. Given the large errors, the predictions from all minima are in agreement with the experimental measurement. We note that the prediction from the global minimum is the one that approaches closest to the experimental measurement, albeit only slightly.

Conclusions
We performed a global fit to the available data in b → cτν τ transitions, adopting an EFT approach with a minimal set of assumptions: 1) NP only enters in the third generation of fermions. 2) There is a sizeable energy gap between NP and the electroweak scale, the EFT operators are SU (2) L ⊗ U (1) Y invariant and the electroweak symmetry breaking is linearly realized. 3) All Wilson coefficients are real (CP is conserved). We have tested the impact of the latter assumption, but did not find an improved description of the data.
In contrast to previous works, we considered the q 2 distributions measured by BaBar and Belle. Moreover, we study the effect of including the recently announced F D * L measurement by the Belle collaboration in the fit. A comparison with earlier analyses, either not including the q 2 distributions, the F D * L measurement, or considering smaller sets of operators, precisely illustrates the benefits of our fit: as described in Section 3, most of the NP solutions found in previous fits are disfavoured once all the information considered in this work is added.
We performed the global fit in different scenarios. As a baseline, we considered the full dataset before the announcement of the F D * L measurement with the subset of operators implied by our assumptions, i.e. with a flavour-universal coefficient C V R . We then performed extensive comparisons to datasets including the recent F D * L measurement, the preliminary Belle measurement of R D ( * ) , and different bounds on B(B c → τν τ ), as well as a second parameter set, allowing for a non-universal C V R .
In the baseline fit, three minima have been obtained, given in Table 2. The global minimum, referred to in the text as Min 1, has an excellent χ 2 ; while none of the fitted Wilson coefficients are required to be non-zero for this minimum, the simplest interpretation of this solution is a global modification of the SM: setting all Wilson coefficients but C V L to zero increases the χ 2 only by ∆χ 2 = 1.4, implying an even better fit. The other two solutions are local minima which numerically exhibit stronger deviations from the SM, with larger contributions of the tensor and scalar operators. While the global minimum is compatible with a SM-like scenario, Min 2 and Min 3 require additional operators. For instance, they could involve scalar LQs with quantum numbers R 2 ∼ (3, 2, 7/6) or S 1 ∼ (3, 1, 1/3).
The measurement of the D * longitudinal polarization fraction F D * L has quite a strong impact on our EFT analysis. It removes Min 3 as a solution for the fit, which was, however, already strongly disfavoured by the differential distributions. Fig. 7 illustrates the tension between the present measurement of F D * L , the bound on B(B c → τν τ ), and the observation ∆R D * > 0 : the set of operators considered within our assumptions cannot accommodate all three observations at 1σ for any combination of Wilson coefficients. Indeed, including the F D * L measurement in the fit increases the minimal χ 2 significantly also for the two lower-lying minima (Min 1b and Min 2b), see Table 3.
We find that most of the minima saturate the upper bound B(B c → τν τ ) ≤ 10%, and it is interesting to study the effect of changing this constraint on the fit. As shown in Tables 2 and 3, adopting a more conservative upper bound of B(B c → τν τ ) ≤ 30% we find the same number of minima; they are qualitatively similar to the previous ones, but with larger central values and ranges of the scalar Wilson coefficients, specifically their pseudoscalar combination. While even this larger upper bound is saturated in most of our fits, the overall decrease in χ 2 is small.
The fact that F D * L cannot be accommodated within 1σ for C V R = 0 could have important consequences, should the present value be confirmed with higher precision. This led us to investigate the scenario with non-zero C V R as a possible resolution of this tension on the theory side. We find that its inclusion helps to reduce the tension among the experimental B → D * data, and it is now possible to satisfy all constraints at 1σ, as illustrated in Fig. 9. The global fit including C V R leads to four different minima, as Fig. 8 shows. Two of these minima have a significantly lower χ 2 than the previous fits, however, they correspond to finetuned solutions where the SM coefficient becomes very small and its effect is substituted by several sizable NP contributions, especially C V R . This scenario seems therefore not to be a satisfactory resolution of the tension.
We have also presented predictions for selected b → cτν τ observables, such as R Λc , R J/ψ or the forward-backward asymmetries and τ polarization in B → D ( * ) τν τ , which have not been included in the fits because either they have not been measured yet or their current experimental values have too large uncertainties. We have studied these observables for the different solutions emerging from our fits, finding that they provide complementary information to the existing data. This is displayed in Figs. 10, 11 and 12. The future measurement of these observables could both establish NP in these modes and allow for a discrimi-nation among the currently favoured scenarios.
We conclude that the anomaly in b → cτν τ transitions remains and can be addressed by NP contributions. Apart from R D ( * ) , also the differential q 2 distributions, F D * L and B(B c → τν τ ) are important to constrain NP, leaving only two viable minima in the global fit. Our general EFT approach does not allow to identify uniquely the potential mediator, since the global minimum can be generated by several combinations of parameters. The generality of our analysis on the other hand allows to use the obtained parameter ranges in more general SMEFT analyses. An improved measurement of F D * L close to its present central value holds the exciting potential to invalidate this general approach, which would have major implications, like a Higgs sector different from the SM one, the existence of NP particles relatively close to the electroweak scale, or new light degrees of freedom. As we have shown, additional measurements will be able to clarify these questions.

D UV Lagrangian
Possible new mediators contributing to the effective Hamiltonian of Eq. (21) and their relative effective Lagrangian are summarized in Table 11.

E Warsaw basis
The operators describing the SMEFT in the Warsaw basis are given by [68,69], O (3) lequ = ¯ j σ µν e ε jk q k σ µν u , where τ I are the Pauli matrices and ε jk is the totally antisymmetric tensor with ε 12 = +1. The fields q and are the quark and lepton SU (2) L doublets, respectively, and u, d, e are the right-handed SU (2) L singlets. Neglecting the small corrections proportional to the CKM factors V ub and V cb , the relevant contributions to the b → cτ ν transitions originate in the Wilson coefficients [C (3) lq ] 3323 ≡C V L , [C (1) lequ ] 3332 ≡C S R , [C ledq ] 3332 ≡C S L and [C (3) lequ ] 3332 ≡C T , where [C X ] ijkl denotes the coefficient of the corresponding operator O X with flavour indices i, j, k, l. The effective Lagrangian relevant for the description of the B anomalies is therefore given by Notice that there is a correspondence between the effective operators at the SMEFT basis with those at the WET basis, according to: which allow us to use the notationC i for the Wilson coefficients at the SMEFT basis, with the aim of making the discussion more intuitive for the reader.
Spin NP mediator Contribution Relevant effective Lagrangian (+ h.c.)  (3), SU (2), U (1) Y ), contribution to the EFT operators and their relevant effective Lagrangian after integrating them out are described for each new field. Their SU (2) decomposition is explicitly shown after the "∼".