On lepton flavour universality in semileptonic $B_c \to \eta_c, J/\psi$ decays

We discuss $B_c \to \eta_c$ and $B_c \to J/\psi$ semileptonic decays within the Standard Model (SM) and beyond. The relevant transition form factors, being the main source of theoretical uncertainties are calculated in the light-cone sum rules approach taking into account all the uncertainties of the model and are provided in a full $q^2$ range. We calculate the semileptonic branching fractions and find for the ratios, $R_{\eta_c}|_{\rm SM}= 0.32 \pm 0.02$, $R_{J/\psi}|_{\rm SM} = 0.23 \pm 0.01$. Both predictions are in agreement with other existing calculations and support the current tension in $R_{J/\psi}$ at 2$\sigma$ level with the experiment. The SM is extended by assuming the most general effective Hamiltonian describing the $b\rightarrow c\tau\nu$ transition, using the constraints on the new physics couplings from the recent combined analysis of BaBar, Belle and LHCb data on semileptonic $B \to D^{(\ast)}$ decays, where the effects of new physics could be visible. We find that these different new physics scenarios are not able to explain the current experimental value of $R_{J/\psi}$. To extend the potential of testing the SM in the semileptonic $B_c$ decays, we consider the forward-backward asymmetry and polarization observables. We also study the 4-fold differential distributions of $B_c \to J/\psi (J/\psi \to \tilde{\ell}^-\tilde{\ell}^+ ) \ell^- \bar{\nu}_\ell$, where $\tilde{\ell} = e, \mu$, in the presence of different new physics scenarios and find that the new physics effects can significantly modify the angular observables and can also produce effects which do not exist in the SM.


Introduction
In September 2017 the LHCb collaboration announced the first measurement of testing the lepton flavor universality using charmed-beauty meson semileptonic decays to J/ψτ + ν µ and J/ψµ + ν µ [1]. The result for the measurement of the ratio of the branching fractions is R J/ψ | exp = BR(B + c → J/ψτ + ν τ ) BR(B + c → J/ψµ + ν µ ) = 0.71 ± 0.17 ± 0.18, (1.1) and is more than 2σ away from the Standard model (SM) prediction. Currently there are many model dependent calculations of R J/ψ [2][3][4][5][6][7][8][9][10][11][12][13][14] within the SM and they give the results in the range (without including model uncertainties) However, the R J/ψ measurement is very challenging. Due to the presence of invisible ν's, both decays are observed only through 3 muons, two of them coming from J/ψ decays and being perfectly identified. The third muon makes a difference and enables distinguishing the semileptonic B c decays to τ and to µ from the background. Therefore it is still premature to speak about the new physics effects in these decays, although one can consider this probability having in mind that BABAR, Belle and LHCb have also found other intriguing anomalies in the semileptonic decays of B mesons, known as R D and R D * [15][16][17][18]. These experimental collaborations have revealed a significant deviation of 2.3σ, and 3.5σ of the ratios R D and R D * from the SM predictions. Also some deviations in b → s semileptonic decays are still present [19]. Moreover, calculations of these semileptonic heavy-meson decays involve theoretical uncertainties coming from imprecise determination of the hadronic transition form factors describing the hadronic effect in the transition from the initial to the final meson state.
The calculation of B c form factors are difficult and leads to big uncertainties. If we summarize values of B c into S-wave charmonia form factors at q 2 = 0 calculated in different models (perturbative QCD (pQCD) [2], three-point QCD sum rules (3ptQCDSR) [3,4], light cone sum rules (LCSR) [5], relativistic quark model (RQM) [6][7][8][9], nonrelativistic quark models (NRQM) [10,11], light-font quark model (LFQM) [12], constituent quark model (CQM)) [13], relativistic quark model (RCQM) [14]) in the literature we obtain: for B c → J/ψ form factors. It is obvious that with such a large range of estimated form factors it is impossible to make any reliable prediction for R ηc and R J/ψ ratios. Moreover, in many estimations of form factors, the theoretical errors were not given or they are not under control. Although some of the uncertainties cancel in the ratio, the model predictions of R J/ψ calculated in different approaches and taking the theoretical uncertainties into account vary in a huge range [20][21][22][23][24][25] R J/ψ | SM = 0.17 − 0.41. (1.5) The lattice QCD calculation for B c → J/ψ form factors V (q 2 ) and A 1 (q 2 ) are available now from the preliminary results of the HPQCD collaboration, at several points for V (q 2 ) and A 1 (q 2 ) [26]. Earlier, the same collaboration has also produced results for B c → η c form factors, which were reported on in the same proceedings. In this paper we will address the calculation of the form factors for B c → S-wave charmonia in the full q 2 range using the LCSR approach. The LCSR method was proven to be a reliable method for calculating transition form factors of many heavy-to light decays, such as B (s) , D (s) → π, ρ, K, K * , η, η [27][28][29][30][31] and even for Λ b → Λ c decays [32,33]. We will compare our results with the existing QCD lattice points for f + (q 2 ) and f 0 (q 2 ) from B c → η c and V (q 2 ) and A 1 (q 2 ) from B c → J/ψ and will show the nice agreement, specially having in mind that the lattice results are still preliminary and do not include systematical errors. Following [22,23] we have assigned 20% uncertainty to the lattice QCD results.
We will also cite the results of recent calculation derived by using the 3ptQCDSR method [4] and will show that the form factors from two sum rule approaches, although calculated by using different quark-hadron duality assumptions, appear to be consistent and precise enough to enable precise determination of the ratios R ηc and R J/ψ in the SM. Recently, there also appeared a modelindependent estimation of the SM bounds on R ηc and R J/ψ [22][23][24]. Such analysis rely on the data, available lattice results and the heavy-quark spin-symmetry (HQSS) relations for the form factors at the zero recoil and predict the R-ratios consistent with Eq.(1.2) and our calculation, clearly at 2σ discrepancy with the experiment. The HQSS and nonrelativistic QCD (NRQCD) relations used in these approaches will be carefully examined in Sec. 2.2.1.
The possible new physics (NP) effects in the semileptonic B c → η c (J/ψ) ν decays have been recently considered in the context of either specific models, such as the leptoquark models [34], left-right symmetric models, R-parity violating supersymmetric models, etc. [35][36][37] or in a model independent approach based on the most general effective Hamiltonian [20,21,25,38]. To account for possible NP effects in B c → η c and B c → J/ψ semileptonic decays we consider here the effective Hamiltonian approach consisting of all possible four-Fermi operators. The constraints on contributions of these NP operators and the corresponding Wilson coefficients are obtained from the experimental results of R D , R D * , polarizations of τ and D * in B → D ( * ) lν decays, as well as on the B c lifetime. There are various studies [39][40][41] performing a global fit on these NP operators considering the presence of only one or two NP operators simultaneously. We have taken the latest constraints on the Wilson coefficients from Ref. [41] and analysed the effects of these NP operators on various observables such as the ratio of the branching fractions, the forward-backward asymmetry, the convexity parameter and the longitudinal as well as the transverse polarization components of τ in the final state. We have also preformed the first study of the full 4-fold differential decay rate B c → J/ψ (J/ψ → µ + µ − , e + e − )lν l , where the leptons from the J/ψ decay are of opposite helicities. The 4-fold decay distribution in this case is proportional to three angles and the momentum transfer q 2 . The three distinct angles gives the freedom to construct additional asymmetries which are sensitive to the real as well as the imaginary part of the new physics couplings.
The structure of the paper is as follows. We compute the form-factors in the context of the LCSR in Sec. 2 and present our results in the whole q 2 range. We compare the results with those existing in the literature and with available preliminary lattice results. The discussion of the heavyquark symmetry limit of form factors at the zero recoil is given in Sec. 2.2.1. The general effective Hamiltonian of the b → c ν transition is introduced in Sec. 3, and we obtain the semileptonic decay distributions for B c → η c , J/ψ in the presence of NP operators using the helicity technique. We compare predictions for different physical observables in the SM and in the presence of NP. A detailed comparison of predictions of R ηc,J/ψ in the SM, with the form factors calculated in LCSR, with the predictions from other approaches is also provided. In Sec. 4, we extend the calculation of B c → J/ψlν l to the J/ψ decay into a pair of muons or electrons, and discuss the full 4-fold distribution. A set of new observables is considered and the results are compared for the SM case and beyond. Finally we conclude in the last section, Sec. 5.

Light-Cone Sum rules calculations for the form factors
We will perform the estimation of the B c → η c and B c → J/ψ form factors using the LCSR method. We will follow the standard LCSR method, by interpolating the B c meson with an appropriate quark current and describing the S-wave charmonia by the light-cone distribution amplitudes (LCDAs) of an increasing twist.
The method of the LCSR is very well know and we will just briefly outline the procedure here in order to properly define all ingredients necessary for calculating the form factors. In the calculation we will use the following approximations: the twist-2 light-cone distribution amplitudes will be calculated in the NRQCD model [42], and the Gegenbauer polynomials expanded at the scale µ. The twist-3 and twist-4 LCDA will be taken in their asymptotic form. Moreover, the Wandzura-Wilczek approximation will be applied, where the three-particle DA are neglected and therefore the twist-3 and twist-4 LCDA are expressed in terms of the twist-2 distributions. The effects of the final state masses, m ηc and m J/ψ are included [43].
The calculation of the form factors for B c → η c proceeds in a similar way as those for B → π, K. [28,30,31,[44][45][46][47], while B c → J/ψ form factor calculation closely follows the derivation of the form factors of B → K * [29,43,46,48,49]. We have checked that with appropriate changes in the expressions, all our results agree with previous calculations.

Definitions
The form factors of the B c → η c decay are defined as where f + (0) = f 0 (0) and 0 ≤ q 2 ≤ (m Bc − m ηc ) 2 . The scalar form factor f 0 (q 2 ) follows also from the conservation of the vector current as The decay B c → J/ψ + ν l is described by the following form factors defined as [29] J/ψ(p, where is the polarization vector of the J/ψ meson, q = p Bc − p is the momentum transfer varying in the range 0 ≤ q 2 ≤ (m Bc − m J/ψ ) 2 and satisfying the relation The form factor A 0 is the pseudoscalar form factor which can also be defined by applying the equation of motion to the derivative of the axial current: and, as it can be seen below, contributes to the B c → J/ψlν decay only if the lepton in the decay is considered to have a non-vanishing mass, which will be a case for the τ particle. The tensor form factors are usually defined as and However, as discussed in [29,43], in the LCSR one has to consider the off-shell p Bc momentum (p 2 Bc = m Bc ) and in order to avoid any ambiguity in the interpretation of p 2 Bc appearing at different steps of calculation it is more appropriate to use the following matrix element as a definition of the tensor form factors; where A(q 2 ), B(q 2 ) and C(q 2 ) are related to T i (q 2 ) defined in Eq. (2.7) as In the LCSR approach the form factors are extracted from the correlation function of the Tproduct of the weak current j Γ , Γ = V, A, S, P, T and an interpolating current of the B c meson j Bc = m bc iγ 5 b among the vacuum and the external on-shell meson M (M = J/ψ, η c ), Both B c → M decays proceed through b-quark decays and we assume that in the region of the large virtualities, the correlation function Eq. (2.11) are dominated by the light-like distances and the description in terms of the products of perturbatively calculable hard-scattering kernels with non-perturbative and universal LCDA, ordered by increasing twist, is appropriate. By inserting the sum over states with B c quantum numbers and by using for the ground state, with the use of hadronic dispersion relation in the virtuality p 2 Bc of the B c channel, we can relate the correlation function Eq. (2.11) to the B c → M matrix elements and the form factors defined above. As usual, the quark-hadron duality is used to approximate heavier state contribution by introducing the effective threshold parameter s Bc 0 and the ground state contribution of the B c meson is enhanced by the Borel transformation in the variable p 2 Bc → σ 2 . The strategy which we use to fix the LCSR parameters, in particular the continuum threshold parameter s Bc 0 is to use the lattice results for the decay constant of B c , Eq. (2.18) and fix the continuum threshold parameters by calculating the constant with the 2-point functions calculated in the LCSR. This is done by using the NLO expression and the pole m b , m c masses. The MS masses used in the paper are taken as m b (m b ) = 4.18 GeV and m c (m c ) = 1.27 GeV. We have achieved the stability of the sum rules, i.e. that continuum and higher-order corrections are suppressed and that also the mass of B c is correctly reproduced for µ = 3.9 ± 0.3 GeV. With the calculated s Bc 0 = 46.8 ± 0.8 GeV 2 we have also checked the stability of the LCSR for B c → M transitions. In both cases, the results are very stable on the variation of the Borel parameter, allowing σ 2 to vary between 70 − 90 GeV 2 with almost no change. Other parameters used in the paper are taken from the lattice results or from the NRQCD models as described below. The leading twist-2 LCDA of a η c meson can be defined as follows [50] 0|c while for the J/ψ we have is a gauge integral. In above ξ = u−(1−u), u is a fraction of a longitudinal momentum of a meson M carried by a c-quark and (1 − u) is a fraction of momentum carried by the c-antiquark. The DAs is defined at a scale µ at which the transverse momenta are integrated up to and all momenta below are included in the nonperturbative DAs φ. Other highertwist amplitudes and higher-order corrections are defined similarly. For all details see, for example [43]. The vector and tensor decay constants f J/ψ and f T J/ψ are defined as 0|c(0)γ µ c(0)|J/ψ(p, e (λ) ) = f J/ψ m J/ψ e (λ) µ , 0|c(0)σ µν c(0)|J/ψ(p, e (λ) ) = if T J/ψ (µ)(e (λ) µ p ν − e (λ) ν p µ ), (2.15) where f T J/ψ is renormalization scale dependent: (2.16) and β 0 = 11 − 2/3n f , n f being the number of flavors involved. The decay constant for η c is defined correspondingly as For the decay constants we will use the lattice results f Bc = 0.427(6)(2) GeV [51,52], f ηc = 0.3947(24) GeV [54], (2.18) while for f T J/ψ we will use the value extracted from the ratio derived by considering combined QCDSR and lattice results [55]. The predictions for charmonia decay constants in [55] also nicely agree with the lattice results above.
The Eq. (2.20) can be inverted to give the coefficients of the conformal expansion a P n (µ) = 2(2n + 3) 3(n + 1)(n + 2) and with the help of these coefficients at some low-energy scale µ 0 , the LCDAs φ P (ξ, µ) can be reconstructed at any scale µ.
The distribution amplitudes can also be defined with the help of calculated moments of LCDA at some scale µ as Charmonia particles are flavor-symmetric and therefore their LCDAs are symmetric around u = 1/2. The second moment is calculated in NRQCD [56,57] where µ 0 ∼ 1GeV. The values of the nonrelativistic speeds v 2 of quarks in η c and J/ψ mesons are obtained in NRQCD by including the first-order α s corrections and non-perturbative contributions proportional to v 2 in the analysis of Γ(η c → γγ) and Γ(J/ψ → e + e − ) rates, respectively, and v 2 J/ψ = 0.225 +0.106 −0.088 [42], v 2 ηc = 0.226 +0.123 −0.098 [42], v 2 J/ψ = v 2 ηc = 0.21 ± 0.02 [57,58] have been extracted. As stated in [59], the two-loop [60,61] and three-loop [62] perturbative corrections to the NRQCD predictions for the Γ(J/ψ → e + e − ) decay rate is known to be large. In [63] and [64] the O(v 2 ) and O(α s ) corrections to twist-2 LCDAs of η c and J/ψ have been calculated. At leading order approximation in relative velocity v there is no difference between η c and J/ψ mesons and the results for the moments obtained are valid for both charmonia LCDA. Based on the power-counting rules of NRQCD one would naively expect that v 2 ∼ 0.3. Taking all above into account, we will use the latest improved value [42,65,66] for both charmonia: For the model of twist-2 LCDA at µ 0 = 1 GeV we adopt the Gaussian model [59]: 28) where N σ ≈ 1 is the normalization constant defined from We also use the Wandzura-Wilczek approximation where three-particle twist-3 LCDAs containing quarks and a gluon are neglected. In that case the twist-3 LCDAs of η c are fixed to their asymptotic forms including mass corrections [67]: For the J/ψ meson the situation is somewhat more complicated. In the Wandzura-Wilczek approximation where three-particle LCDAs are neglected, by using equations of motion the twist-3 LCDAs can be expressed in terms of the leading twist-2 DAs φ ,⊥ with the valence quark mass corrections as [27,[68][69][70][71]:h (s) ⊥ , and The J/ψ twist-4 LCDA will be taken in their asymptotic form: Some comments are in order. In the B c → η c decay we will retain only contributions up to twist-3 terms. It is well known that the standard twist expansion works very well for B → pseudoscalar form factors. It could be that in our case, for B c → η c , twist-4 corrections are somewhat larger, due to the large and non-negligible m c mass, but since hadronic parameters for the twist-4 contribution for η c are not known we will not include them. In the decay B c → J/ψ we keep all contributions up to twist-4, since their asymptotic form does not depend on the hadronic parameters. The J/ψ DAs defined above do not correspond to matrix elements of operators with definite twist [29]: φ ⊥, are of twist-2, h (s,t) and g (v,a) ⊥ contain a mixture of twist-2 and 3 contributions and A ⊥, , h 3 , g 3 are a mixture of twist-2, 3 and 4 contributions. Therefore it is usual to refer to g , h (s,t) as twist-3 LCDAs and to h 3 , g 3 , A ⊥, as twist-4 LCDAs. Also, as the mass of the vector particle in B → vector decays plays a significant role, in [29] it was proposed the following classification of relevance in the two-particle LCDA: δ ∼ m J/ψ is treated as an expansion parameter. For more detailed discussion see [43,46,48].

Parametrization of the form factors and the results
The derivation of the LCSR expressions for the form factors proceeds in a standard way [43,46,48,49]. The general expression for the calculation of the form factors in the LCSR approach is given by σ is a Borel parameter andū = 1 − u. The functions F (u, µ, q 2 ) contain all twist contributions in terms of the various twist LCDAs, and the higher-twist contributions are suppressed by the Borel parameter.
The derived results at q 2 = 0 are listed in Table 1, together with the earlier results found in the literature on the same form factors. The errors are obtained by varying all parameters in a given range and adding them in quadratures. Note that we completely disagree with a previous LCSR calculation of the form factors in [5].
[26]  Table 1: Form factor predictions at q 2 = 0. Recent relevant lattice results are given by the HPQCD collaboration [26], reported here in orange, without the systematical error.
It is well know that the form factors extracted from the LCSR are valid in the low q 2 region. We use our LCSR results for the form factors and calculate them in the range q 2 = {−5, 5} GeV. Then we extrapolate them from the low q 2 region to the q 2 max by using Bourrely-Caprini-Lellouch (BCL) parametrization [73] of the form factor series expansion in powers of a conformal mapping variable, which satisfies unitarity, analyticity and perturbative QCD scaling. The BCL parametrization is based on a rapidly converging series in the parameter z as weighted by a simple pole function P (q 2 ) = 1 − t/m 2 R which accounts for low-laying resonances present below the threshold production of real is a free parameter that can be used to optimize the convergence of the series expansion. For the truncation to only two terms in expansion Eq. (2.41), it was shown that the optimized value of t 0 has the form [74]: and that the other choices of t 0 do not make a visible change in the form factors parametrization.
Masses of resonances appearing in the fits are determined by the properties of the form factors. The form factors V and T 1 correspond to the vector components of the currents, and, as the B c meson is a pseudoscalar, they correspond to the axial vector components of the matrix elements. A 1,2 , as well as T 2,3 , correspond to the axial vector component of the V −A, while the form factor A 0 correspond to the pseudoscalar current and only contributes in the decays with the non-vanishing lepton masses (in the semileptonic B c decays with the τ lepton in our case). All relevant resonance masses are given in Table 2, together with the fitted parameters α 0 , α 1 from Eq. (2.41).
The predicted form factors in a full q 2 range are shown in Figs. 1, 2, 3.   Figure 1: Pseudoscalar form factors for Bc → ηc calculated in this paper, including the lattice points from [26] with added 20% systematical error.

HQSS/NRQCD symmetry relations among from factors at the zero recoil
It is interesting to check the HQSS and NRQCD limits of the form factors for B c → η c and B c → J/ψ decays. These decays are specific since the decay proceeds through b → c quark decay and the produced final state is a particle formed by two c quarks. And although they look like as a heavy-to-heavy transitions and produce interesting symmetries in a heavy-quark limit [79], the c-quark is significantly lighter than b and the produced c-quark is quite energetic, which spoils exact heavy-flavor symmetries. On the other hand c-quark is heavy enough that such decays can be considered as nonrelativistic so that the approximation of the zero-recoil point, i.e. the symmetry relations for a maximum momentum transfer q 2 max = (m Bc − m J/ψ,ηc ) 2 still hold and the form factors can be related to a single function ∆ [79][80][81], with an unknown normalization. Following [79] we write for the form factors near zero recoil (q m c ):  Figure 2: SM form factors for Bc → J/ψ calculated in this paper, including the lattice points from [26] with added 20% systematical error. where V µ =bγ µ c and A µ =bγ µ γ 5 c and µ is a polarization vector of J/ψ. Here v is the velocity of the B c meson, and q is a small residual velocity carried by the final state meson (not to be confused by q, the momentum carried by the lepton pair system), so that The parameter a 0 is connected to the Bohr radius of the B c meson, it value is not important for the further discussion and will not be discussed here.
We can now relate the ∆(a 0 q ) function to the B c → η c form factor f + (q 2 ) at the zero recoil as which amounts, using the predicted f + (q 2 max ) from the LCSR calculation above, to ∆(a 0 q ) LCSR ≈ 0.79 ± 0.09. (2.47) This value can be compared with the value obtained in the QCD relativistic potential model in [81].
In [80] it was shown that in the NRQCD approximation one can derive a generalized set of relations using the HQSS, so that the transition form factors of B c → η c and B c → J/ψ decays can be given in terms of a single form factor, even for the case of non-equal four-velocities v 1 = v 2 , of the initial and the final state heavy mesons. If the following helicity basis for the form factors in B c → J/ψ decay is defined the expressions from [80], stemming from considering NRQCD and HQSS, and relating different decay form factors at the point of zero recoil q 2 max ofQq →Q q transitions can be expressed as [22]: for the B c → J/ψ decay, and for the B c → η c decay [23], where some shorthand notation has been introduced: Additionally, the vector decay form factors can be related to the pseudoscalar ones as [24] f (q 2 max ) = where we have used that r q = r Q for B c → η c , J/ψ and simplified the relations. It is expected that these relations are broken by terms of order O(m c /m b , Λ QCD /m c ) 30%.
Here we check the consistency of these relations and our LCSR results. At the zero-recoil F 1 (q 2 ) and f (q 2 ) are the same up to a constant factor, Eq. (2.50), which explicitly gives and We see that the HQSS/NRQCD predictions are quite consistent with the exact LCSR predictions for the form factors at the zero recoil and can be safely used in model-independent bounds on R ratios as it was done in [22][23][24], keeping in mind that their accuracy is limited to O(30%). Finally, using Eq. (2.52) we obtain, 3 R ηc , R J/ψ and decay distributions of B c → η c ν and B c → J/ψ ν The general effective Lagrangian for the quark level transition b → c ν with = e, µ, τ is given by with the four-Fermi operators defined as We use σ µν = i [γ µ , γ ν ] /2 and V L,R , S L,R , T L are the complex Wilson coefficients governing the NP contributions which are zero in the SM. Since we want to explain the possible lepton-flavour non-universality, we will assume that the NP effects contribute to the τ leptons only. The matrix element of the semileptonic decays B c → J/ψ(η c )τ ν τ then has the form: We note that the axial and the pseudoscalar hadronic currents do not contribute to the B c → η c decay, and therefore V R − V L = 0, S R − S L = 0, ⇒ V R = V L , S R = S L . The scalar hadronic current does not contribute to the B c → J/ψ transition which leads to S L + S R = 0. We henceforth use the shorthand definition S R + S L = S and S R − S L = P in the text. The constraints on the Wilson coefficients appearing in Eq. (3.1) are obtained from the combined analysis of the BaBar, Belle and LHCb data for the branching fraction ratios R D ( * ) , the τ polarization asymmetry along the longitudinal directions of the τ lepton in B → D * , as well as the longitudinal D * polarization in B c → D * τ ν τ decay [41]. The leptonic branching fraction of the B c meson, BR(B c → τ ν), is not yet measured, therefore the possible NP contributions come from precise experimental measurements of the B c lifetime, τ exp Bc = (0.507 ± 0.009) ps [82]. The theoretical SM prediction of the B c lifetime still allows for up to 60% contribution from NP [41,83] in the B c leptonic decay width. In particular, the best fit point for S R is dependent on the assumption of the B c → τ ν decay width.
We consider for our analysis the limit BR(B c → τν) < 30% and the values of the Wilson coefficients from the combined analysis done in Ref. [41]. They studied all one-dimensional scenarios with only one NP Wilson coefficient considered at a time and the two-dimensional scenarios with two NP Wilson coefficients considered simultaneously. The best fit points in the 1D scenarios and their 2σ ranges (given in square brackets below) at 1 TeV are given in Table 1  All NP operators are generated by the addition of a single new particle to the SM. The relation S L = 4 T L is generated in the R 2 leptoquark scenario with a scalar SU (2) L doublet [85,86] at the new physics scale. The leptoquark model with an SU (2) L singlet scalar S 1 gives the relation S L = −4 T L at the NP scale. These relations are modified at the scale m b to: , after including one-loop electroweak corrections in addition to the three-loop QCD anomalous dimensions in the renormalization group running using the following relations [87]: We now discuss the differential decay rates for the processes B c → η c ν and B c → J/ψ ν . The differential decay rate for these semi-leptonic processes depend on the angle θ which is the polar angle of the lepton (the angle between the lepton direction in the W * rest frame and the direction of the W * in the B c rest frame) and the momentum transfer q 2 ( q = p Bc − p) to the ν pair. The differential (q 2 , cos θ ) distribution can be calculated using the helicity techniques and is of the form where |p 2 | = λ 1/2 (m 2 Bc , m 2 ηc,J/ψ , q 2 )/2m Bc is the momentum of η c (J/ψ) in the B c rest frame, v = (1−m 2 /q 2 ) is the lepton velocity in the −ν center-of-mass frame and H µν L µν is the contraction of the hadronic and the leptonic tensors. The helicity techniques to calculate the angular distribution in the presence of new physics operators for the semi-leptonic decays considered here can be found in Ref. [88,89].
The differential distribution for the B c → η c τ ν τ is written as with P = p Bc + p and q = p Bc − p (p = p ηc or p = p J/ψ . Next, the differential distribution of theB c → J/ψ −ν decay is considered with V R V L = ReV R ReV L + ImV R ImV L , T L P = ReT L ReP + ImT L ImP and P V L = ReP ReV L + ImP ImV L , and is given by The hadronic helicity amplitudes in terms of the form factors given in Eqs. (2.3, 2.7) are expressed as 3.1 Results for the branching ratios and R ηc , R J/ψ predictions We first give our predictions for branching fractions in the SM of both decays in Table 3 Table 3: Branching fractions of Bc → J/ψ, ηc decays calculated in different models and given in %, with l denoting a light lepton, e or µ The ratios of semileptonic branching fractions using our form factors derived by the LCSR calculation expressed in Eq. (2.39) gives We see that above results agree with the recent model-independent analysis of R J/ψ [22,24] and R ηc [23,24]. See also the discussion in Sec.2.2.1.
The ratios of the branching fractions R J/ψ,ηc are computed next in the context of different NP scenarios using the form factors calculated using LCSR in Sec. 2. The values of the NP operators' effective couplings considered for our analysis are discussed before and are given by Eqs. (3.4, 3.5). In Fig. 4 we show the the q 2 dependence of the ratios R ηc and R J/ψ in the presence of only one NP operator (first two figures of both panels). The third figure in both panels shows the ratio in presence of two NP operators. The SM value is always shown by the blue dotted line. We see that the ratio increases for most of NP contributions for both J/ψ and η c . The S L = 4 T L case with the coupling being pure real or imaginary results in a decrease in the ratio R ηc . This is due to the negative interference between S L and T L , Eq. (3.8). The shaded region shows the 2σ allowed region for V L , S L = 4T L , S L,R parameters in the 1D fit, with the central value shown by a dashed line. In case of the 2D scenarios the results are presented for the best fit point. As expected, the ratio R ηc is more sensitive to the scalar and the tensor operators, whereas R J/ψ is more sensitive to V L . The values of R J/ψ and R ηc in the presence of different NP scenarios are listed in Table 4. The results are presented for the best fit points, as well as for the 2σ allowed regions in the 1D scenario. Note that any of the considered NP scenarios derived from the recent global fit analysis on available experimental data on semileptonic B → (D, D * ) ν decays [41] cannot explain the 2σ tension with the experiment Eq. 1.1 of R J/ψ ratio.

Forward-backward asymmetry, convexity parameter and the τ polarization
The differential distributions defined in Eqs. (3.8, 3.10) can be written in a simple form as a function of cos θ as

Bc
(A(q 2 ) + B(q 2 ) cos θ + C(q 2 ) cos 2 θ ). (3.14) Observables depending on the polar angle distribution of the emitted leptons such as the forwardbackward lepton asymmetry and the convexity parameter are considered first. They are defined by The A(q 2 ), B(q 2 ) and C(q 2 ) functions can be easily obtained from Eqs. (3.8, 3.10). We present in Figs. 5, 6 the q 2 dependence of the forward-backward asymmetry A F B (q 2 ) and the convexity parameter C τ F (q 2 ). These observables are not sensitive to the case where V L is the only NP contribution. The B c → η c transition appears to be more sensitive to the new physics operators as compared to the B c → J/ψ transition. In case of the J/ψ decay mode, the presence of the S L , S R coefficients in the 2D scenario leads to a significant deviation from A F B (q 2 ) prediction in the SM. The present allowed values of the coupling have a very small effect on C τ F (q 2 ) in case of J/ψ, whereas in case of η c the S L = 4 T L case enhances C τ F (q 2 ) only at large values of q 2 .  Now we discuss the effect on the polarization of the emitted τ in the W − rest frame in the presence of the NP operators. The differential decay rate for a given spin projection in a given direction can be easily obtained with the inclusion of the spin projection operators (1 + γ 5 / s i )/2 for τ in the calculation. The longitudinal and the transverse polarization components of the τ are then defined as:

1D Scenario
where s µ L and s µ T are the longitudinal and the transverse polarization four-vectors of τ − in the W − rest frame and are given by [90][91][92] s µ L = 1 m τ (| p τ |, E τ sin θ τ , 0, E τ cos θ τ ), s µ T = (0, cos θ τ , 0, − sin θ τ ).  The longitudinal and transverse polarizations in the B c → η, J/ψτ ν τ decays are given as : The transverse polarization of τ as can be seen from Eq. (3.19) has an overall factor of √ δ τ and therefore vanishes in the limit of zero lepton mass and the emitted lepton is then fully longitudinally polarized. Therefore, the τ lepton can be largely transversely polarized as compared to the muons or the electrons. The q 2 dependence of the τ polarization in presence of different NP operators is shown in Figs. 7, 8. The following observations can be made from the figures. The longitudinal and transverse polarizations of τ in the η c decay mode are more sensitive to the NP operators compared to the J/ψ decay mode. The tau transverse polarization in the J/ψ decay mode is again mostly affected by the NP operator S L = 4 T L at low values of q 2 , whereas the S L , S R parameters in the 2D  scenario lead to a deviation from the SM prediction for both the longitudinal and the transverse τ polarization. The predictions for the mean forward-backward asymmetry, the convexity parameter and the tau polarization in the presence of different NP operators are summarised in Table 5.  The blue dotted lines are the SM prediction, the green dashed line is for the best fit values of the NP couplings in the 1D scenario as discussed in the text. The green band represents the NP effects from the 2σ allowed regions. The third figure in both panels is the result for the best fit points in the 2D scenarios. 4 Decay distribution of B c → J/ψ (J/ψ → µ + µ − ) ν decay We consider in this section the process B c → J/ψ (J/ψ → µ + µ − ) ν , with the 4-fold differential decay rate being dependent on three angles θ V , θ , χ and the momentum transfer q 2 . The angle θ is same as defined before, θ V is the polar angle between the direction of the emitted µ − in the J/ψ rest frame and the parent J/ψ in the B c rest frame, and χ is the azimuthal angle between the W * ν plane and the J/ψµ + µ − plane. The angles are shown in Fig. 9 and are defined as usually being taken in the literature. The J/ψ is too light to decay to τ + τ − , therefore the outgoing leptons can be either a pair of µ or of e. We ignore the mass m µ , m e from the J/ψ decays but the mass of lepton from W * decay is retained. The total differential decay rate for the µ − L µ + R (σ ∼ λ − − λ + = −1) final state is given by Eq. (4.1) below. The corresponding expressions for µ − R µ + L final state can be obtained by setting θ V → θ V + π in Eq. 4.1.
ReP cos χ + ImP ) sin χ − H −− cos 2 θ V 2 ReP cos χ − ImP sin χ , We only list the interference terms of NP with SM and do not show the NP-NP interference terms, but they are included in our calculations. The expressions above are now more involved for the 4-fold differential distribution and contain various combinations with θ , θ V and χ angles, with the imaginary couplings being proportional to sin χ. The constraints on the NP coefficients (V L , S L and S R ) in the 1D scenario are obtained using the condition that they are purely real. The global fit results which we consider here, do not include the vector operator V R , as this vector operator with right-handed coupling to the quarks does not arise at the dimension-six level in the SU (2) Linvariant effective theory. The relation S L = 4 T L in the pure imaginary case is in more agreement with the SM compared to the case with the real Wilson coefficients. However, the effects of the real and the imaginary components of these NP coefficients can be isolated by constructing different angular asymmetries. We first consider the forward-backward asymmetry in θ V and both θ V , θ with the angle χ fully integrated over: where G = d 4 Γ dq 2 d cos θ d cos θ V dχ and Γ in the denominator is the decay width of B c → µ + µ − ν , obtained by integrating Eq. (4.1) and is given by It can be seen from Eq. (4.2) that the numerator is not sensitive to the scalar type NP operators. Therefore the sensitivity to the scalar NP comes only from the total decay width in the denominator, Eq. (4.3). In Fig. 10  One can build additional asymmetries in the angle χ along with θ V and θ . These asymmetries   are proportional to both cos χ and sin χ, and their corresponding expressions are given below: We show in Fig. 11, the asymmetries A Finally, we consider the observables which will be sensitive to only the imaginary component of the NP operators. These asymmetries are zero within the SM. There are three possible combinations, (a) asymmetry depending only on χ, (b) asymmetry depending on χ and θ V , and (c)   asymmetry depending on χ, θ V and θ . The relevant expressions are: where (T L V L ) * = ImT L ReV L −ImV L ReT L , (P V L ) * = ImP ReV L −ImV L ReP, (P T L ) * = ImP ReT L − ImT L ReP , The asymmetry A img F B (χ, θ V , θ ) in only sensitive to the NP operator V R and is therefore not relevant for our case since these V R coefficients are not considered in the global fits as discussed before. We show in Fig. 12  These observables are only shown for ImS L = 4 ImT L in the 1D scenario and ReS L = 4 ReT L , ImS L = 4 ImT L in the 2D scenarios as these were the only cases considered in the global fit in Ref. [41]. The forward-backward asymmetry depending only on χ in the light of results from the current global fit shows about 1% deviation from the SM in the 1D scenario and up to 3% deviation in the 2D scenario, in the mid-range of q 2 = 5 − 9 GeV 2 . The asymmetry A img F B (χ, θ V ) in case of the 2D scenario will have around only 1% deviation in low q 2 region, whereas it is insensitive to NP in the 1D scenario. The predictions for the integrated forward-backward asymmetries in the presence of  Table 6: The integrated values of the forward-backward asymmetries in the whole q 2 region, in case of different NP scenarios discussed in the text. The subscript and the superscript are the values for the 2σ range of the NP couplings.

Conclusions
Experimental measurements of semileptonic decays of the B mesons have lead to intriguing experimental tensions with the SM in the last years. The LHCb measurement of B c → J/ψlν l decays has lead to the speculation whether the observed potential lepton flavour universality (LFU) violation in B decays can be also seen in the semileptonic B c channels. However, the SM prediction for the B c decays require a knowledge of the transition form factors of B c → η c , J/ψ and the ignorance of the form factor theoretical errors yields a degree of uncertainty in the prediction. Preliminary results for these form factors exist at couple of q 2 values from the lattice QCD, but they do not cover the entire allowed range of the momentum transfer and are still given without systematical errors. We have calculated the form factors in the LCSR approach and have given the results in the full q 2 region. Our results are in good agreement with the existing lattice points. The SM branching ratios of the B c meson to J/ψ and η c are calculated in LCSR and compared with the results from other approaches. Our predictions for the semileptonic ratios R J/ψ | SM = 0.23 ± 0.01 and R ηc | SM = 0.32 ± 0.02 are in agreement with other derivations and support the existing tension at 2σ level with the experiment on R J/ψ , Eq. (1.1). With more data on B c decays from HL-LHC all observables in B c → η c , J/ψ semileptonic decays will be within the reach of LHCb and tested in the near future.
The possible NP effects in the semileptonic decays of B c to η c and J/ψ is also studied based on the effective Hamiltonian approach consisting of all possible four-fermi operators. The constraints on these NP operators can be obtained from the experimental data on R D( * ) , the τ and D * longitudinal polarization from B → D * decay and the leptonic B c → τ µ branching ratio. We take into account the latest constraints from Ref. [41] and analyse the effects of the NP operators on various observables. The ratio R J/ψ is sensitive to V L in the high q 2 range whereas R ηc is more sensitive to the scalar and the tensor operators, as expected.
The sensitivity of all the considered observables in this work to the different NP operators is summarized in Table 7. We find that most of the observables in the η c decay mode are sensitive to the NP coupling S R . The transverse polarization of τ is mostly affected by the current best fit point of the combination of coefficients Re,Im[S L = 4 T L ] in the 2D NP scenario. The 2D NP scenario with the presence of both S R , S L leads to the largest deviation from the SM predictions for most of the observables in the case of J/ψ, apart from R J/ψ .  Table 7: Summary of the sensitivity of the observables to the NP couplings. The best fit value of the NP coupling which is most sensitive to the observable is marked with . The boxes with * are the ones where 2σ ranges of NP parameters give the largest deviation from the SM value.
In addition, the full 4-fold differential distribution of the decay rate B c → J/ψ ν , with J/ψ decaying to a pair of leptons of opposite helicity is considered for the fist time in the presence of new physics operators. We find that the asymmetry in the angle θ V (A J/ψ F B (θ V )) is mostly sensitive to the NP couplings Re,Im[S L = 4T L ], in the 2D NP scenarios. The asymmetries in the angle χ, which are zero in the SM and are sensitive to the imaginary part of the NP coupling, are also considered and found to be sensitive to S L = 4 T L combination of parameters. Therefore, with the current allowed parameter space for the S L = 4 T L NP parameters obtained from the global fit to experimental data on semileptonic B → D, D * decays, the asymmetries constructed with θ V , χ and (θ V , χ) angles lead to significant deviation from the SM prediction.
However, it is important to stress that none of the NP scenarios derived from the recent global fit analysis of the available experimental data on semileptonic B → (D, D * ) ν decays [41] can also simultaneously explain the current 2σ tension with the experimental R J/ψ ratio. With the extended experimental LHCb program, future studies with more data will boost or disapprove this evidence of LFU violation in B c decays.