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

TL;DR

This paper analyzes semileptonic $B_c$ decays to $ a_c$ and $J/\psi$ within the Standard Model and beyond, calculating form factors, branching ratios, and exploring new physics effects on angular observables, highlighting current experimental tensions.

Contribution

It provides detailed form factor calculations over the full $q^2$ range and investigates new physics impacts on angular distributions in $B_c$ decays, extending the testing framework beyond existing studies.

Findings

Predicted $R_{ a_c}|_{SM} = 0.32 \\pm 0.02$ and $R_{J/\psi}|_{SM} = 0.23 \\pm 0.01.

Current new physics scenarios cannot fully explain the experimental $R_{J/\psi}$ value.

New physics can significantly alter angular observables and produce effects absent in the SM.

Abstract

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 sum rule approach 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. 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 \rightarrow 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. 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 simultaneously explain the current experimental value of $R_{J/\psi}$.

Tables

Table 1. Table 1 : Form factor predictions at q 2 = 0 superscript 𝑞 2 0 q^{2}=0 . Recent relevant lattice results are given by the HPQCD collaboration [ 26 ] , reported here in orange, without the systematical error.

Form Factor	this work 222The cited errors are obtained just by varying all parameters of the model and their smallness shows the stability of the sum rules used to obtain the predictions for the form factors. The errors do not include intrinsic uncertainties of the model itself which are hard to predict and could potentially increase the errors and lower the accuracy of the predictions.	QCDSR [4]	QCDSR  [3]	SR [5]	pQCD  [2]	CCQM  [25, 77]	RQM  [9]	RQM  [8]	LFQM  [12]	latt. [26]
$f_{+,0}^{{\eta }_c}\left(0\right)$	$0.62\pm 0.05$	$0.41\pm 0.04$	0.66	0.87	0.48(7)	0.75	0.47	0.54	0.61(5)	0.59
$V^{J/\psi }\left(0\right)$	$0.73\pm 0.06$	$0.70\pm 0.06$	1.03	1.69	0.42(2)	0.78	0.49	0.73	0.74(4)	0.70
$A_1^{J/\psi }\left(0\right)$	$0.55\pm 0.04$	$0.50\pm 0.05$	0.63	0.75	0.46(3)	0.56	0.50	0.52	0.50(3)	0.48
$A_2^{J/\psi }\left(0\right)$	$0.35\pm 0.03$	$0.43\pm 0.05$	0.69	1.69	0.64(3)	0.55	0.73	0.51	0.44(5)	-
$A_0^{J/\psi }\left(0\right)$	$0.54\pm 0.04$	$0.53\pm 0.04$	0.60	0.27	0.59(3)	0.56	0.40	0.53	0.53(3)	-
$f_T^{{\eta }_c}\left(0\right)$	$0.93\pm 0.07$	-	-	-	-	0.93	-	-	-	-
$T_{1,2}^{J/\psi }\left(0\right)$	$0.47\pm 0.04$	$0.48\pm 0.03$	-	-	-	0.56	-	-	-	-
$T_3^{J/\psi }\left(0\right)$	$0.19\pm 0.01$	$0.27\pm 0.03$	-	-	-	0.20	-	-	-	-

Table 2. Table 2 : Summary of the BCL fit for B c → η c → subscript 𝐵 𝑐 subscript 𝜂 𝑐 B_{c}\to\eta_{c} and B c → J / ψ → subscript 𝐵 𝑐 𝐽 𝜓 B_{c}\to J/\psi form factors. The masses of the low-laying B c subscript 𝐵 𝑐 B_{c} resonances are taken from [ 80 , 81 , 82 , 83 ]

Form factor	$J^P$	$m_R$ (GeV)	${\alpha }_0$	${\alpha }_1$
$f_{+}$	$1^{-}$	6.34	0.62	-6.13
$f_0$	$0^{+}$	6.71	0.63	-4.86
$f_T$	$1^{-}$	6.34	0.93	-9.36
$V$	$1^{-}$	6.34	0.74	-8.66
$A_1$	$1^{+}$	6.75	0.55	-4.67
$A_2$	$1^{+}$	6.75	0.35	-1.78
$A_0$	$0^{-}$	6.28	0.54	-6.80
$T_1$	$1^{-}$	6.34	0.48	-4.88
$T_2$	$1^{+}$	6.75	0.48	-2.93
$T_3$	$1^{+}$	6.75	0.19	-1.69

Table 3. Table 3 : Branching fractions of B c → J / ψ , η c → subscript 𝐵 𝑐 𝐽 𝜓 subscript 𝜂 𝑐 B_{c}\rightarrow J/\psi,\eta_{c} decays calculated in different models and given in %, with l 𝑙 l denoting a light lepton, e 𝑒 e or μ 𝜇 \mu .

Mode	this work	QCDSR  [3]	SR [5]	pQCD  [2]	RCQM  [14]	CCQM  [25, 77]	RQM  [9]	RQM  [6]	RQM  [8]	RQM  [7]	LFQM  [12]
$B_c\to {\eta }_{c}l{\overline{\nu }}_l$	${0.82}_{-0.11}^{+0.12}$	0.85 (0.75)	1.85 (1.64)	0.50 (0.44)	0.91 (0.81)	0.95	0.47 (0.42)	0.89	0.52 (0.52)	0.85	0.74 (0.67)
$B_c\to {\eta }_c\tau {\overline{\nu }}_{\tau }$	${0.26}_{-0.05}^{+0.06}$	0.25 (0.23)	0.55 (0.49)	0.15 (0.14)	0.25 (0.22)	0.24	-	-	-	-	0.21 (0.19)
$B_c\to J/\psi l{\overline{\nu }}_l$	${2.24}_{-0.49}^{+0.57}$	2.16 (1.9)	2.67 (2.37)	1.13 (1.00)	2.33 (2.07)	1.67	1.39 (1.23)	1.42	1.49 (1.47)	2.33	1.64 (1.49)
$B_c\to J/\psi \tau {\overline{\nu }}_{\tau }$	${0.53}_{-0.14}^{+0.16}$	0.54 (0.48)	0.73 (0.65)	0.33 (0.29)	0.55 (0.49)	0.40	-	-	-	-	0.41 (0.37)

Equations

$$R_{J/\psi }{|}_{\exp }=\frac{BR(B_c^{+}\to J/\psi {\tau }^{+}{\nu }_{\tau })}{BR(B_c^{+}\to J/\psi {\mu }^{+}{\nu }_{\mu })}=0.71\pm 0.17\pm 0.18,$$

$$R_{J/\psi }{|}_{\mathrm{SM}}=0.24-0.30,$$

$$\displaystyle f_{+}(0)=f_{0}(0)=0.20-1.43,\$$

$$V(0)=0.17-1.63,$$

$$A_1(0)=0.21-1.19,$$

$$A_2(0)=0.23-1.27,$$

$$A_0(0)=0.12-1.09,$$

$$R_{J/\psi }{|}_{\mathrm{SM}}=0.17-\mathrm{0.41.}$$

$$⟨{\eta }_c(p)|\bar{c}{\gamma }_{\mu }b|B_c(p_{B_c})⟩$$

$$⟨{\eta }_c(p)|\bar{c}{\sigma }_{\mu \nu }{q}^{\nu }b|B_c(p_{B_c})$$

$$⟨0|\bar{c}b|B_c(p_{B_c})⟩=\frac{m_{B_c}^2-m_{{\eta }_c}^2}{m_b(\mu )-m_c(\mu )}{f}_0(q^2).$$

$$⟨J/\psi (p,ϵ)|\bar{c}{\gamma }_{\mu }(1-{\gamma }_5)b|B_c(p_{B_c})⟩=-i{ϵ}_{\mu }^{\ast }(m_{B_c}+m_{J/\psi })A_1(q^2)+i(p_{B_c}+p{)}_{\mu }({ϵ}^{\ast }\cdot q)\frac{A_2(q^2)}{m_{B_c}+m_{J/\psi }}$$

$$A_3(q^2)$$

$$A_3(0)=A_0(0).$$

$$⟨J/\psi |\bar{c}i{\gamma }_{5}b|B_c⟩=\frac{2m_{J/\psi }}{(m_b(\mu )+m_c(\mu ))}({ϵ}^{\ast }\cdot q)A_0(q^2)$$

$$⟨J/\psi (p,ϵ)|\bar{c}{\sigma }_{\mu \nu }{q}^{\nu }(1+{\gamma }_5)b|B(p_{B_c})⟩=2i{\epsilon }_{\mu \nu \rho \sigma }{ϵ}^{\ast \nu }{p}_{B_c}^{\rho }{p}^{\sigma }{T}_1(q^2)$$

$$T_1(0)=T_2(0).$$

$$⟨J/\psi (p,ϵ)|\bar{c}{\sigma }_{\mu \nu }{\gamma }_{5}b|B(p_{B_c})⟩$$

$$T_1(q^2)=A(q^2),T_2(q^2)=A(q^2)-\frac{q^2}{m_{B_c}^2-m_{J/\psi }^2}B(q^2),T_3(q^2)=B(q^2)+C(q^2).$$

$$\Pi (q^2,p_{B_c}^2)$$

$$⟨0|j_{B_c}|B_c(p_B)⟩=\frac{f_{B_c}{m}_{B_c}^2}{m_b(\mu )+m_c(\mu )}$$

$$⟨0|\bar{c}(z){\gamma }_{\mu }{\gamma }_5[z,-z]c(-z)|{\eta }_c(p)⟩=-i{f}_{{\eta }_c}{p}_{\mu }{\int }_{-1}^{1}d\xi e^{i\xi pz}\varphi (\xi ,\mu ),$$

$$⟨0|\bar{c}(z){\gamma }_{\mu }[z,-z]c(-z)|J/\psi (p,{ϵ}^{(\lambda =0)})⟩=f_{J/\psi }{m}_{J/\psi }{p}_{\mu }{\int }_{-1}^{1}d\xi e^{i\xi px}{\varphi }_{||}(\xi ,\mu ),$$

$$⟨0|\bar{c}(z){\sigma }_{\mu \nu }[z,-z]c(-z)|J/\psi (p,{ϵ}^{(\lambda =\pm 1)})⟩=i{f}_{J/\psi }^{\perp }({ϵ}_{\mu }{p}_{\nu }-{ϵ}_{\nu }{p}_{\mu }){\int }_{-1}^{1}du{e}^{i\xi px}{\varphi }_{\perp }(\xi ,\mu ),$$

$$⟨0|\bar{c}(0){\gamma }_{\mu }c(0)|J/\psi (p,e^{(\lambda )})⟩$$

$$⟨0|\bar{c}(0){\sigma }_{\mu \nu }c(0)|J/\psi (p,e^{(\lambda )})⟩$$

$$f_{J/\psi }^T({\mu }^{\prime 2})={\left({\alpha }_s({\mu }^{\prime 2})/{\alpha }_s({\mu }^2)\right)}^{C_f/{\beta }_0}{f}_{J/\psi }^T(\mu )$$

$$⟨0|\bar{c}{\gamma }_{\mu }{\gamma }_{5}c|{\eta }_c(p)⟩=-i{f}_{{\eta }_c}{p}_{\mu }.$$

$$f_{B_c}$$

$$f_{J/\psi }$$

$$f_{{\eta }_c}$$

$$R_{J/\psi }^T=\frac{f_{J/\psi }^T(\mu =2\mathrm{GeV})}{f_{J/\psi }}=0.975\pm 0.010,$$

$${\varphi }_P(u,{\mu }^2)=6u(1-u)\left(1+\sum _{n=1}^{\infty }{a}_n^P({\mu }^2)C_n^{3/2}(2u-1)\right).$$

$$a_n^P(Q^2)={\left({\alpha }_s(Q^2)/{\alpha }_s({\mu }^2)\right)}^{{\gamma }_n^P/(2{\beta }_0)}{a}_n^P({\mu }^2),$$

Full text

1 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/\psi\tau^{+}\nu_{\mu}$ and $J/\psi\mu^{+}\nu_{\mu}$ [1]. The result for the measurement of the ratio of the branching fractions is

[TABLE]

and is more than 2$\sigma$ away from the Standard model (SM) prediction. Currently there are many model dependent calculations of $R_{J/\psi}$ [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)

[TABLE]

However, the $R_{J/\psi}$ measurement is challenging. Due to the presence of invisible $\nu$’s, both decays are observed only through 3 muons, two of them coming from $J/\psi$ decays and being perfectly identified. The third muon makes a difference and enables distinguishing the semileptonic $B_{c}$ decays to $\tau$ and to $\mu$ 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^{\ast}}$ [15, 16, 17, 18]. These experimental collaborations have revealed a significant deviation of 2.3$\sigma$, and 3.5$\sigma$ of the ratios $R_{D}$ and $R_{D^{\ast}}$ from the SM predictions. Also some deviations in $b\to 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:

[TABLE]

for $B_{c}\to\eta_{c}$ form factors, and

[TABLE]

for $B_{c}\to J/\psi$ 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_{\eta_{c}}$ and $R_{J/\psi}$ 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/\psi}$ calculated in different approaches and taking the theoretical uncertainties into account vary in a huge range [20, 21, 22, 23, 24, 25]

[TABLE]

The lattice QCD calculation for $B_{c}\to J/\psi$ 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}\to\eta_{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}\to S$-wave charmonia in the full $q^{2}$ range using the LCSR-inspired 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)}\to\pi,\rho,K,K^{\ast},\eta,\eta^{\prime}$ [27, 28, 29, 30, 31] and even for $\Lambda_{b}\to\Lambda_{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}\to\eta_{c}$ and $V(q^{2})$ and $A_{1}(q^{2})$ from $B_{c}\to J/\psi$ 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 also cite the results of our recent calculation derived by using the 3ptQCDSR method [4]111A very brief discussion on the 3ptQCDSR calculation is provided in Appendix A. and 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_{\eta_{c}}$ and $R_{J/\psi}$ in the SM. Recently, there also appeared a model-independent estimation of the SM bounds on $R_{\eta_{c}}$ and $R_{J/\psi}$ [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.(2) and our calculation, clearly at $2\sigma$ 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}\to\eta_{c}(J/\psi)\ell\nu_{\ell}$ 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, 39, 40]. To account for possible NP effects in $B_{c}\to\eta_{c}$ and $B_{c}\to J/\psi$ 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^{\ast}}$, polarizations of $\tau$ and $D^{\ast}$ in $B\to D^{(\ast)}l\nu$ decays, as well as on the $B_{c}$ lifetime. There are various studies [41, 42, 43, 44, 45, 46] 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. [43] 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 $\tau$ in the final state. We have also preformed the first study of the full 4-fold differential decay rate $B_{c}\to J/\psi\,(J/\psi\to\mu^{+}\mu^{-},e^{+}e^{-})l\nu_{l}$, where the leptons from the $J/\psi$ 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 our sum rule model in Sec. 2 and present the 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 heavy-quark symmetry limit of form factors at the zero recoil is given in Sec.  2.2.1. The general effective Hamiltonian of the $b\rightarrow c\ell\nu_{\ell}$ transition is introduced in Sec. 3, and we obtain the semileptonic decay distributions for $B_{c}\to\eta_{c},J/\psi$ 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_{\eta_{c},J/\psi}$ in the SM, with the form factors calculated in our model, with the predictions from other approaches is also provided. In Sec. 4, we extend the calculation of $B_{c}\to J/\psi l\nu_{l}$ to the $J/\psi$ 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.

2 Sum rule calculations for the form factors

We will perform the estimation of the $B_{c}\to\eta_{c}$ and $B_{c}\to J/\psi$ form factors using the LCSR-inspired method. We will follow the standard QCD sum rule method, by interpolating the $B_{c}$ meson with an appropriate quark current and describing the S-wave charmonia by the distribution amplitudes (DAs) 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 [47], and the Gegenbauer polynomials expanded at the scale $\mu$. The twist-3 and twist-4 DAs 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 DAs are expressed in terms of the twist-2 distributions. The effects of the final state masses, $m_{\eta_{c}}$ and $m_{J/\psi}$ are included [48].

The calculation of the form factors for $B_{c}\to\eta_{c}$ proceeds in a similar way as those for $B\to\pi,K.$ [28, 30, 31, 49, 50, 51, 52], while $B_{c}\to J/\psi$ form factor calculation closely follows the derivation of the form factors of $B\to K^{*}$ [53, 51, 29, 54, 48]. We have checked that with appropriate changes in the expressions, all our results agree with previous calculations.

2.1 Definitions

The form factors of the $B_{c}\to\eta_{c}$ decay are defined as

[TABLE]

where $f_{+}(0)=f_{0}(0)$ and $0\leq q^{2}\leq(m_{B_{c}}-m_{\eta_{c}})^{2}$. The scalar form factor $f_{0}(q^{2})$ follows also from the conservation of the vector current as

[TABLE]

The decay $B_{c}\to J/\psi\ell^{+}\nu_{l}$ is described by the following form factors defined as [29]

[TABLE]

where $\epsilon$ is the polarization vector of the $J/\psi$ meson, $q=p_{B_{c}}-p$ is the momentum transfer varying in the range $0\leq q^{2}\leq(m_{B_{c}}-m_{J/\psi})^{2}$ and

[TABLE]

satisfying the relation

[TABLE]

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:

[TABLE]

and, as it can be seen below, contributes to the $B_{c}\to J/\psi l\nu$ decay only if the lepton in the decay is considered to have a non-vanishing mass, which will be a case for the $\tau$ particle. The tensor form factors are usually defined as

[TABLE]

and

[TABLE]

However, as discussed in [29, 48], in the standard QCD sum rule one has to consider the off-shell $p_{B_{c}}$ momentum ($p_{B_{c}}^{2}\neq m_{B_{c}}$) and in order to avoid any ambiguity in the interpretation of $p_{B_{c}}^{2}$ appearing at different steps of calculation it is more appropriate to use the following matrix element as a definition of the tensor form factors;

[TABLE]

where $A(q^{2})$, $B(q^{2})$ and $C(q^{2})$ are related to $T_{i}(q^{2})$ defined in Eq. (2.1) as

[TABLE]

The form factors are extracted from the correlation function of the T-product of the weak current $j_{\Gamma},\Gamma=V,A,S,P,T$ and an interpolating current of the $B_{c}$ meson $j_{B_{c}}=m_{b}\bar{c}i\gamma_{5}b$ among the vacuum and the external on-shell meson M ($M=J/\psi,\eta_{c}$),

[TABLE]

Both $B_{c}\to M$ decays proceed through $b$-quark decays and we assume that in the region of the large $m_{b}^{2}-q^{2}\leq\mathcal{O}(m_{b}\Lambda_{QCD})$ and $m_{b}^{2}-p_{B_{c}}^{2}\leq\mathcal{O}(m_{b}\Lambda_{QCD})$ virtualities, the correlation function Eq. (16) 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 light-cone distribution amplitude (LCDA), ordered by increasing twist, is appropriate.

By inserting the sum over states with $B_{c}$ quantum numbers and by using

[TABLE]

for the ground state, with the use of hadronic dispersion relation in the virtuality $p_{B_{c}}^{2}$ of the $B_{c}$ channel, we can relate the correlation function Eq. (16) to the $B_{c}\to 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_{0}^{B_{c}}$ and the ground state contribution of the $B_{c}$ meson is enhanced by the Borel transformation in the variable $p_{B_{c}}^{2}\to\sigma^{2}$.

The strategy which we use to fix the sum rule parameters, in particular the continuum threshold parameter $s_{0}^{B_{c}}$, is to use the lattice results for the decay constant of $B_{c}$, Eq. (2.1) 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 $\overline{\rm MS}$ masses used in the paper are taken as $m_{b}(\overline{m_{b}})=4.18$ GeV and $m_{c}(\overline{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 $\mu=3.9\pm 0.3$ GeV. With the calculated $s_{0}^{B_{c}}=46.8\pm 0.8$ GeV2 we have also checked the stability of the sum rules for $B_{c}\to M$ transitions. In both cases, the results are very stable on the variation of the Borel parameter, allowing $\sigma^{2}$ to vary between $70-90\,{\rm GeV}^{2}$ with almost no change. Other parameters used in the paper are taken from the lattice results or from the NRQCD models described afterwards.

The method of the LCSR was extensively used for calculating the heavy-to-light transition form factors. Here the situation is far more complicated since the final meson is a quarkonium state $\eta_{c}$ or $J/\psi$. Therefore, to properly account for the non-negligible large mass corrections ${\cal O}(2m_{c})$ in the correlator, one would have to do a systematic expansion of the correlator near the light cone including those corrections. This is a highly nontrivial task and according to our knowledge has not been done yet. In the future, to improve the whole picture, one has to do a revised consideration of the LCDA for charmonia, similar to what was done for heavy hadrons ($B$-mesons and $\Lambda_{b}$), by proving the factorization theorems and deriving the RG evolution kernels of LCDA by considering full mass corrections. But, such a calculation for charmonia is far more complicated since there is no help from HQET and heavy-quark symmetries, neither can one achieve fast convergence in the heavy-mass expansion. Such a calculation, if consistently doable for charmonia, is far beyond the scope of our paper. Here we assume that these potentially large intrinsic mass effects can be effectively described using proper phenomenological model of DAs. So, we will follow a simplified sum rule model where we treat charmonia of $B_{c}$-decays as light particles in the correlator (16) expansion near the light-cone and will closely follow the approach of the standard LCSR in what follows. On the other hand, to describe nonperturbative properties of charmonia we will use the NRQCD-inspired DAs which exactly reproduce leading NR moments of charmonia at $\sim$ 1 GeV energies. To resolve the right DA structure at the $\sim m_{b}$ energies of the decay, we calculate the Gegenbauer expansion and evolution of DAs. The corrections to the leading approximation will be done by making the twist-expansion and by taking the large $m_{c}$ mass correction in twist-3 and twist-4 DAs into account. The genuine $\mathcal{O}(4m_{c}^{2}/m_{b}\Lambda_{QCD})$ corrections are not included as we assume the collinearity of the wave functions. Moreover, since we are aware of our model constraints in describing charmonia particles, we will show the stability of the model on the variation of parameters of the model, the consistency of our results with the 3ptSR calculation of the same form factors done with the same parameters used here and will also show consistency of the calculated form factors with the HQSS/NRQCD symmetry relations among form factors.

The leading twist-2 DA of a $\eta_{c}$ meson is defined as follows [55]

[TABLE]

while for the $J/\psi$ we have

[TABLE]

where $[z,-z]=P\exp{ig\int_{-z}^{z}dx^{\mu}A_{\mu}(x)}$ is a gauge integral. In above $\xi=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 $\mu$ at which the transverse momenta are integrated up to and all momenta below are included in the nonperturbative DAs $\phi$. Other higher-twist amplitudes and higher-order corrections are defined similarly. For all details see, for example [48]. The vector and tensor decay constants $f_{J/\psi}$ and $f_{J/\psi}^{T}$ are defined as

[TABLE]

where $f_{J/\psi}^{T}$ is renormalization scale dependent:

[TABLE]

and $\beta_{0}=11-2/3n_{f}$, $n_{f}$ being the number of flavors involved. The decay constant for $\eta_{c}$ is defined correspondingly as

[TABLE]

For the decay constants we will use the lattice results

[TABLE]

while for $f_{J/\psi}^{T}$ we will use the value extracted from the ratio

[TABLE]

derived by considering combined QCDSR and lattice results [60]. The predictions for charmonia decay constants in [60] also nicely agree with the lattice results above.

2.1.1 Distribution amplitudes for charmonia

The leading twist-2 DAs are expanded in terms of Gegenbauer polynomials as:

[TABLE]

The leading term is the asymptotic form $\phi(u,\mu^{2}\to\infty)=6u(1-u)$. The Gegenbauer coefficients $a_{n}$ are renormalized multiplicatively

[TABLE]

where the anomalous dimensions $\gamma^{\parallel,\perp}_{n}$ are given by

[TABLE]

Here $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ and $\beta_{0}=11/3N_{c}-2/3N_{f}$, in which $N_{c}$ is the number of colors and $N_{f}$ the number of flavors. The coefficients $a_{n}^{||}$ appear in $\phi(\xi,\mu)$ and $\phi_{||}(\xi,\mu)$, while $a_{n}^{\perp}$ are coefficients in the expansion of the transversal twist-2 DA $\phi_{\perp}(\xi,\mu)$ of a $J/\psi$ meson.

The Eq. (25) can be inverted to give the coefficients of the conformal expansion

[TABLE]

and with the help of these coefficients at some low-energy scale $\mu_{0}$, the DA $\phi_{P}(\xi,\mu)$ can be reconstructed at any scale $\mu$.

The distribution amplitudes can also be defined with the help of calculated moments of DAs at some scale $\mu$ as

[TABLE]

Charmonia particles are flavor-symmetric and therefore their DAs are symmetric around $u=1/2$. The second moment is calculated in NRQCD [61, 62]

[TABLE]

where $\mu_{0}\sim 1{\rm GeV}$.

The values of the nonrelativistic speeds $v^{2}$ of quarks in $\eta_{c}$ and $J/\psi$ mesons are obtained in NRQCD by including the first-order $\alpha_{s}$ corrections and non-perturbative contributions proportional to $v^{2}$ in the analysis of $\Gamma(\eta_{c}\to\gamma\gamma)$ and $\Gamma(J/\psi\to e^{+}e^{-})$ rates, respectively, and $\langle v^{2}\rangle_{J/\psi}=0.225\,_{-0.088}^{+0.106}$ [47], $\langle v^{2}\rangle_{\eta_{c}}=0.226\,_{-0.098}^{+0.123}$ [47], $\langle v^{2}\rangle_{J/\psi}=\langle v^{2}\rangle_{\eta_{c}}=0.21\pm 0.02$ [63, 62] have been extracted. As stated in [64], the two-loop [65, 66] and three-loop [67] perturbative corrections to the NRQCD predictions for the $\Gamma(J/\psi\to e^{+}e^{-})$ decay rate is known to be large. In [68] and [69] the ${\cal O}(v^{2})$ and ${\cal O}(\alpha_{s})$ corrections to twist-2 DAs of $\eta_{c}$ and $J/\psi$ have been calculated. At leading order approximation in relative velocity $v$ there is no difference between $\eta_{c}$ and $J/\psi$ mesons and the results for the moments obtained are valid for both charmonia DAs. Based on the power-counting rules of NRQCD one would naively expect that $\langle v^{2}\rangle\sim 0.3$. Taking all above into account, we will use the latest improved value [47, 70, 71] for both charmonia:

[TABLE]

For the model of twist-2 DA at $\mu_{0}=1$ GeV we adopt the Gaussian model [64]:

[TABLE]

where $N_{\sigma}\approx 1$ is the normalization constant defined from

[TABLE]

We also use the Wandzura-Wilczek approximation where three-particle twist-3 DAs containing quarks and a gluon are neglected. In that case the twist-3 DAs of $\eta_{c}$ are fixed to their asymptotic forms including mass corrections [72]:

[TABLE]

where $\rho_{+}(\mu)=4m_{c}^{2}(\mu)/m_{\eta_{c}}^{2}$ and

[TABLE]

For the $J/\psi$ meson the situation is somewhat more complicated. In the Wandzura-Wilczek approximation where three-particle DAs are neglected, by using equations of motion the twist-3 DAs can be expressed in terms of the leading twist-2 DAs $\phi_{\parallel,\perp}$ with the valence quark mass corrections

[TABLE]

as [27, 73, 74, 75, 76]:

[TABLE]

and

[TABLE]

with

[TABLE]

The $J/\psi$ twist-4 DAs will be taken in their asymptotic form:

[TABLE]

Some comments are in order. In the $B_{c}\to\eta_{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\to$ pseudoscalar form factors. It could be that in our case, for $B_{c}\to\eta_{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 $\eta_{c}$ are not known we will not include them. In the decay $B_{c}\to J/\psi$ we keep all contributions up to twist-4, since their asymptotic form does not depend on the hadronic parameters. The $J/\psi$ DAs defined above do not correspond to matrix elements of operators with definite twist [29]: $\phi_{\perp,\parallel}$ are of twist-2, $h_{\parallel}^{(s,t)}$ and $g_{\perp}^{(v,a)}$ contain a mixture of twist-2 and 3 contributions and ${\mathbb{A}}_{\perp,\parallel}$, $h_{3}$, $g_{3}$ are a mixture of twist-2, 3 and 4 contributions. Therefore it is usual to refer to $g_{\perp}^{(v,a)},h_{\parallel}^{(s,t)}$ as twist-3 LCDAs and to $h_{3},g_{3},{\mathbb{A}}_{\perp,\parallel}$ as twist-4 LCDAs. Also, as the mass of the vector particle in $B\to$ vector decays plays a significant role, in [29] it was proposed the following classification of relevance in the two-particle LCDA: $O(\delta^{0}):$ $\phi_{\perp}$; $O(\delta^{1}):$ $\phi_{\parallel},g_{\perp}^{(v,a)}$; $O(\delta^{2}):$ $h_{\parallel}^{(s,t)},h_{3},{\mathbb{A}}_{\perp}$; $O(\delta^{3}):$ $g_{3},{\mathbb{A}}_{\parallel}$, where now $\delta\sim m_{J/\psi}$ is treated as an expansion parameter. For more detailed discussion see [51, 53, 48].

2.2 Parametrization of the form factors and the results

The derivation of the sum rule expressions for the form factors proceeds in a standard way [53, 51, 54, 48].

The general expression for the calculation of the form factors is given by

[TABLE]

where

[TABLE]

$\sigma$ is a Borel parameter and $\bar{u}=1-u$. The functions $F(u,\mu,q^{2})$ contain all twist contributions in terms of the various twist DAs, 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 recent QCDSR result [4] briefly discussed in Appendix A and 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.

It is well know that the form factors extracted from the sum rules are valid in the low $q^{2}$ region. We use our 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_{\rm max}^{2}$ by using Bourrely-Caprini-Lellouch (BCL) parametrization [78] 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

[TABLE]

weighted by a simple pole function $P(q^{2})=1-t/m_{R}^{2}$ which accounts for low-laying resonances present below the threshold production of real $B_{c}-M$ pairs at $t_{+}=(m_{B_{c}}+m_{M})^{2}$. The parameter $t_{0}$, $0\leq t_{0}\leq t_{-}=(m_{B_{c}}-m_{M})^{2}$ 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. (46), it was shown that the optimized value of $t_{0}$ has the form [79]:

[TABLE]

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 $\tau$ lepton in our case). All relevant resonance masses are given in Table 2, together with the fitted parameters $\alpha_{0},\alpha_{1}$ from Eq. (46).

The predicted form factors in a full $q^{2}$ range are shown in Figs. 1, 2, 3.

2.2.1 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}\to\eta_{c}$ and $B_{c}\to J/\psi$ decays. These decays are specific since the decay proceeds through $b\to 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 [84], 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}_{\rm max}=(m_{B_{c}}-m_{J/\psi,\eta_{c}})^{2}$ still hold and the form factors can be related to a single function $\Delta$ [84, 85, 86], with an unknown normalization. Following [84] we write for the form factors near zero recoil ($q^{\prime}\ll m_{c}$):

[TABLE]

where $V_{\mu}=\bar{b}\gamma_{\mu}c$ and $A_{\mu}=\bar{b}\gamma_{\mu}\gamma_{5}c$ and $\epsilon_{\mu}$ is a polarization vector of $J/\psi$. Here $v$ is the velocity of the $B_{c}$ meson, and $q^{\prime}$ 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

[TABLE]

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 $\Delta(a_{0}q^{\prime})$ function to the $B_{c}\to\eta_{c}$ form factor $f_{+}(q^{2})$ at the zero recoil as

[TABLE]

which amounts, using the predicted $f_{+}(q_{\rm max}^{2})$ from the calculation above, to

[TABLE]

This value can be compared with the value obtained in the QCD relativistic potential model in  [86].

In [85] 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}\to\eta_{c}$ and $B_{c}\to J/\psi$ decays can be given in terms of a single form factor, even for the case of non-equal four-velocities $v_{1}\neq v_{2}$, of the initial and the final state heavy mesons. If the following helicity basis for the form factors in $B_{c}\to J/\psi$ decay is defined

[TABLE]

and

[TABLE]

the expressions from [85], stemming from considering NRQCD and HQSS, and relating different decay form factors at the point of zero recoil $q^{2}_{\text{max}}$ of $\bar{Q}q\rightarrow\bar{Q}^{{}^{\prime}}q$ transitions can be expressed as [22]:

[TABLE]

for the $B_{c}\rightarrow J/\psi$ decay, and

[TABLE]

for the $B_{c}\rightarrow\eta_{c}$ decay [23], where some shorthand notation has been introduced: $r_{M}=m_{M}/m_{B_{c}}$ (with $m_{M}=[m_{J/\psi},m_{\eta_{c}}]$), $r_{Q}=m_{Q^{{}^{\prime}}}/m_{Q}=m_{c}/m_{b}$ and $r_{q}=m_{q}/m_{Q}=r_{Q}$. Additionally, the vector decay form factors can be related to the pseudoscalar ones as [24]

[TABLE]

where we have used that $r_{q}=r_{Q}$ for $B_{c}\to\eta_{c},J/\psi$ and simplified the relations. It is expected that these relations are broken by terms of order $\mathcal{O}(m_{c}/m_{b},\Lambda_{\text{QCD}}/m_{c})\lesssim 30\%$.

Here we check the consistency of these relations and our form factor results. At the zero-recoil $\mathcal{F}_{1}(q^{2})$ and $f(q^{2})$ are the same up to a constant factor, Eq. (55), which explicitly gives

[TABLE]

whereas

[TABLE]

and

[TABLE]

We see that the HQSS/NRQCD predictions are quite consistent with our sum rule 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 ${\cal O}(30\%)$. Finally, using Eq. (57) we obtain,

[TABLE]

an excellent agreement among the relations between $B_{c}\to\eta_{c}$ and $B_{c}\to J/\psi$ transition form factors derived from HQSS/NRQCD symmetry relations and our exact results at the zero recoil.

3 $R_{\eta_{c}}$, $R_{J/\psi}$ and decay distributions of $B_{c}\rightarrow\eta_{c}\ell\nu_{\ell}$ and $B_{c}\rightarrow J/\psi\ell\nu_{\ell}$

The general effective Lagrangian for the quark level transition $b\rightarrow c\ell\nu_{\ell}$ with $\ell=e,\mu,\tau$ is given by

[TABLE]

with the four-Fermi operators defined as

[TABLE]

We use $\sigma_{\mu\nu}=i\left[\gamma_{\mu},\gamma_{\nu}\right]/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 $\tau$ leptons only. The matrix element of the semileptonic decays $B_{c}\rightarrow J/\psi(\eta_{c})\tau\nu_{\tau}$ then has the form:

[TABLE]

We note that the axial and the pseudoscalar hadronic currents do not contribute to the $B_{c}\rightarrow\eta_{c}$ decay, and therefore $V_{R}-V_{L}=0,S_{R}-S_{L}=0,\Rightarrow V_{R}=V_{L},S_{R}=S_{L}$. The scalar hadronic current does not contribute to the $B_{c}\rightarrow J/\psi$ 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. (62) are obtained from the combined analysis of the BaBar, Belle and LHCb data for the branching fraction ratios $R_{D^{(\ast)}}$, the $\tau$ polarization asymmetry along the longitudinal directions of the $\tau$ lepton in $B\to D^{\ast},$ as well as the longitudinal $D^{\ast}$ polarization in $B_{c}\rightarrow D^{\ast}\tau\nu_{\tau}$ decay [43]. The leptonic branching fraction of the $B_{c}$ meson, $BR(B_{c}\to\tau\nu)$, is not yet measured, therefore the possible NP contributions come from precise experimental measurements of the $B_{c}$ lifetime, $\tau^{\mathrm{exp}}_{B_{c}}=(0.507\pm 0.009)$ ps [87]. The theoretical SM prediction of the $B_{c}$ lifetime still allows for up to 60% contribution from NP [88, 43] 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}\to\tau\nu$ decay width.

We consider for our analysis the limit BR$(B_{c}\rightarrow\tau\bar{\nu})<$ 30% and the values of the Wilson coefficients from the combined analysis done in Ref. [43]. 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\sigma$ ranges (given in square brackets below) at 1 TeV are given in Table 1 of Ref. [43] and we list them below for completion:

[TABLE]

Only the real values of the coefficients were considered for the fit. The possibility of allowing imaginary coefficients was examined in Ref. [89] and they obtained that the relation ${\rm Im}[S_{L}]=4\,{\rm Im}[T_{L}]$ is also permitted by the recent experiments. We therefore use the best fit value for $S_{L}=4\,T_{L}$ in Eq. (3) for both, the real and the imaginary case. The results of the fit for the NP Wilson coefficients in the 2D scenario at 1 TeV are taken from Table 2 of Ref. [43]:

[TABLE]

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 [90, 91] 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: $S_{L}(m_{b})\simeq 8.1\,T_{L}(m_{b})$ for $R_{2}$, and $S_{L}(m_{b})\simeq-8.5\,T_{L}(m_{b})$, after including one-loop electroweak corrections in addition to the three-loop QCD anomalous dimensions in the renormalization group running using the following relations [92]:

[TABLE]

We now discuss the differential decay rates for the processes $B_{c}\rightarrow\eta_{c}\ell\nu_{\ell}$ and $B_{c}\rightarrow J/\psi\ell\nu_{\ell}$. The differential decay rate for these semi-leptonic processes depend on the angle $\theta_{\ell}$ which is the polar angle of the lepton $\ell$ (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_{B_{c}}-p$) to the $\ell\nu_{\ell}$ pair. The differential $(q^{2},\cos\theta_{\ell})$ distribution can be calculated using the helicity techniques and is of the form

[TABLE]

where $|{\bf p_{2}}|=\lambda^{1/2}(m_{B_{c}}^{2},m_{\eta_{c},J/\psi}^{2},q^{2})/2m_{B_{c}}$ is the momentum of $\eta_{c}(J/\psi)$ in the $B_{c}$ rest frame, $v=(1-m_{\ell}^{2}/q^{2})$ is the lepton velocity in the $\ell^{-}\bar{\nu}_{\ell}$ center-of-mass frame and $H_{\mu\nu}L^{\mu\nu}$ 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. [93, 94].

The differential distribution for the $B_{c}\rightarrow\eta_{c}\tau\nu_{\tau}$ is written as

[TABLE]

with the helicity flip-factor $\delta_{\ell}=m_{\ell}^{2}/2q^{2}$, $T_{L}V_{L}={\rm Re}T_{L}~{}{\rm Re}V_{L}+{\rm Im}T_{L}~{}{\rm Im}V_{L}$, $T_{L}S={\rm Re}T_{L}~{}{\rm Re}S+{\rm Im}T_{L}~{}{\rm Im}S$ and $S\,V_{L}={\rm Re}S~{}{\rm Re}V_{L}+{\rm Im}S~{}{\rm Im}V_{L}$. We consider the interference between the different NP operators as their effect can be similar to NP2, if they are of the same value. The $H^{\prime}s$ in Eq. (69) are the hadronic helicity amplitudes written in terms of the invariant form factors defined in Eq. (2.1) and are of the form

[TABLE]

with $P=p_{B_{c}}+p$ and $q=p_{B_{c}}-p$ ($p=p_{\eta_{c}}$ or $p=p_{J/\psi}$.

Next, the differential distribution of the $\bar{B}_{c}\to J/\psi\ell^{-}\bar{\nu}_{\ell}$ decay is considered with $V_{R}V_{L}={\rm Re}V_{R}~{}{\rm Re}V_{L}+{\rm Im}V_{R}~{}{\rm Im}V_{L}$, $T_{L}P={\rm Re}T_{L}~{}{\rm Re}P+{\rm Im}T_{L}~{}{\rm Im}P$and $P\,V_{L}={\rm Re}P~{}{\rm Re}V_{L}+{\rm Im}P~{}{\rm Im}V_{L}$, and is given by

[TABLE]

The hadronic helicity amplitudes in terms of the form factors given in Eqs. (8, 2.1) are expressed as

[TABLE]

3.1 Results for the branching ratios and $R_{\eta_{c}},R_{J/\psi}$ predictions

We first give our predictions for branching fractions in the SM of both decays in Table 3, where we put the branching fraction values updated using the latest value for the $B_{c}$ lifetime, $\tau_{B_{c}}=(0.507\pm 0.009)$ ps  [87], while in the brackets we cite the original published values of the $BR$s. If there are no brackets the branching fractions have already been calculated using the latest value for $\tau_{B_{c}}$.

The ratios of semileptonic branching fractions using our calculated form factors from Eq. (44) are

[TABLE]

We see that above results agree with the recent model-independent analysis of $R_{J/\psi}$ [22, 24] and $R_{\eta_{c}}$ [24, 23]. See also the discussion in Sec.2.2.1.

The ratios of the branching fractions $R_{J/\psi,\eta_{c}}$ are computed next in the context of different NP scenarios using the form factors calculated in Sec. 2. The values of the NP operators’ effective couplings considered for our analysis are discussed before and are given by Eqs. (3, 3). In Fig. 4 we show the the $q^{2}$ dependence of the ratios $R_{\eta_{c}}$ and $R_{J/\psi}$ 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/\psi$ and $\eta_{c}$. The $S_{L}=4\,T_{L}$ case with the coupling being pure real or imaginary results in a decrease in the ratio $R_{\eta_{c}}$. This is due to the negative interference between $S_{L}$ and $T_{L}$ , Eq. (69). The shaded region shows the $2\sigma$ 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_{\eta_{c}}$ is more sensitive to the scalar and the tensor operators, whereas $R_{J/\psi}$ is more sensitive to $V_{L}$.

The values of $R_{J/\psi}$ and $R_{\eta_{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$\sigma$ 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\to(D,D^{\ast})\ell\nu_{\ell}$ decays [43] cannot explain the 2$\sigma$ tension with the experiment Eq. 1 of $R_{J/\psi}$ ratio.

[TABLE]

3.2 Forward-backward asymmetry, convexity parameter and the $\tau$ polarization

The differential distributions defined in Eqs. (69, 71) can be written in a simple form as a function of $\cos\theta_{\ell}$ as

[TABLE]

Observables depending on the polar angle distribution of the emitted leptons such as the forward-backward lepton asymmetry and the convexity parameter are considered first. They are defined by

[TABLE]

The $\mathcal{A}(q^{2}),\mathcal{B}(q^{2})$ and $\mathcal{C}(q^{2})$ functions can be easily obtained from Eqs. (69, 71). We present in Figs. 5, 6 the $q^{2}$ dependence of the forward-backward asymmetry $A_{FB}(q^{2})$ and the convexity parameter $C^{\tau}_{F}(q^{2})$. These observables are not sensitive to the case where $V_{L}$ is the only NP contribution. The $B_{c}\rightarrow\eta_{c}$ transition appears to be more sensitive to the new physics operators as compared to the $B_{c}\to J/\psi$ transition. In case of the $J/\psi$ decay mode, the presence of the $S_{L},S_{R}$ coefficients in the 2D scenario leads to a significant deviation from $A^{FB}(q^{2})$ prediction in the SM. The present allowed values of the coupling have a very small effect on $C^{\tau}_{F}(q^{2})$ in case of $J/\psi$, whereas in case of $\eta_{c}$ the $S_{L}=4\,T_{L}$ case enhances $C^{\tau}_{F}(q^{2})$ only at large values of $q^{2}$.

Now we discuss the effect on the polarization of the emitted $\tau$ 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+\gamma_{5}\not{s}_{i})/2$ for $\tau$ in the calculation. The longitudinal and the transverse polarization components of the $\tau$ are then defined as:

[TABLE]

where $s^{\mu}_{L}$ and $s^{\mu}_{T}$ are the longitudinal and the transverse polarization four-vectors of $\tau^{-}$ in the $W^{-}$ rest frame and are given by [95, 96, 97]

[TABLE]

The longitudinal and transverse polarizations in the $B_{c}\to\eta,J/\psi\tau\nu_{\tau}$ decays are given as :

[TABLE]

[TABLE]

The transverse polarization of $\tau$ as can be seen from Eq. (80) has an overall factor of $\sqrt{\delta_{\tau}}$ and therefore vanishes in the limit of zero lepton mass and the emitted lepton is then fully longitudinally polarized. Therefore, the $\tau$ lepton can be largely transversely polarized as compared to the muons or the electrons.

The $q^{2}$ dependence of the $\tau$ 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 $\tau$ in the $\eta_{c}$ decay mode are more sensitive to the NP operators compared to the $J/\psi$ decay mode. The tau transverse polarization in the $J/\psi$ 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 $\tau$ 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.

[TABLE]

4 Decay distribution of $B_{c}\rightarrow J/\psi\,(J/\psi\to\mu^{+}\mu^{-})\,\ell\nu_{\ell}$ decay

We consider in this section the process $B_{c}\rightarrow J/\psi\,(J/\psi\rightarrow\mu^{+}\mu^{-})\,\ell\nu_{\ell}$, with the 4-fold differential decay rate being dependent on three angles $\theta_{V},\theta_{\ell},\chi$ and the momentum transfer $q^{2}$. The angle $\theta_{\ell}$ is same as defined before, $\theta_{V}$ is the polar angle between the direction of the emitted $\mu^{-}$ in the $J/\psi$ rest frame and the parent $J/\psi$ in the $B_{c}$ rest frame, and $\chi$ is the azimuthal angle between the $W^{*}\ell\nu$ plane and the $J/\psi\mu^{+}\mu^{-}$ plane. The angles are shown in Fig. 9 and are defined as usually being taken in the literature.

The $J/\psi$ is too light to decay to $\tau^{+}\tau^{-}$, therefore the outgoing leptons can be either a pair of $\mu$ or of $e$. We ignore the mass $m_{\mu},m_{e}$ from the $J/\psi$ decays but the mass of lepton from $W^{*}$ decay is retained. The total differential decay rate for the $\mu^{-}_{L}\mu^{+}_{R}$ ($\sigma\sim\lambda_{\ell}^{-}-\lambda_{\ell}^{+}=-1$) final state is given by Eq. (81) below. The corresponding expressions for $\mu^{-}_{R}\mu^{+}_{L}$ final state can be obtained by setting $\theta_{V}\rightarrow\theta_{V}+\pi$ in Eq. 81.

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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 $\theta_{\ell}$, $\theta_{V}$ and $\chi$ angles, with the imaginary couplings being proportional to $\sin\chi$. 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)_{L}$-invariant 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 $\theta_{V}$ and both $\theta_{V},\theta_{\ell}$ with the angle $\chi$ fully integrated over:

[TABLE]

where $\mathcal{G}=\left(\frac{d^{4}\Gamma}{dq^{2}\,d\cos\theta_{\ell}\,d\cos\theta_{V}\,d\chi}\right)$ and $\Gamma$ in the denominator is the decay width of $B_{c}\to\mu^{+}\mu^{-}\ell\nu_{\ell}$, obtained by integrating Eq. (81) and is given by

[TABLE]

It can be seen from Eq. (4) 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. (83). In Fig. 10 we show $A_{FB}^{J/\psi}(\theta_{V})$ and $A_{FB}^{J/\psi}(\theta_{V},\theta_{\ell})$ as a function of $q^{2}$ with the values of the new physics couplings as given in Eqs. (3, 3). The current bound on the NP couplings makes the observable $A_{FB}^{J/\psi}(\theta_{V})$ sensitive to $S_{L}=4\,T_{L}$ in the 1D scenario and to the same combination with both the real and the imaginary components present in case of 2D scenario. As for the asymmetry $A_{FB}^{J/\psi}(\theta_{V})$, the deviation from the SM in case of 2D scenario for the combination Re[$S_{L}=4\,T_{L}$], Im[$S_{L}=4\,T_{L}$] can be as large as 50-70% in the region of small $q^{2}$. However, the other observable with the asymmetry in both $\theta_{V}$ and $\theta_{\ell}$, $A_{FB}^{J/\psi}(\theta_{V},\theta_{\ell})$ is not a good observable to look for NP scenarios in the current situation.

One can build additional asymmetries in the angle $\chi$ along with $\theta_{V}$ and $\theta_{\ell}$. These asymmetries are proportional to both $\cos\chi$ and $\sin\chi$, and their corresponding expressions are given below:

[TABLE]

We show in Fig. 11, the asymmetries $A_{FB}^{J/\psi}(\chi,\theta_{V})$ [upper-panel] and $A^{FB}(\chi,\theta_{V},\theta_{\ell})$ [lower-panel] as a function of $q^{2}$. They behave similar to the asymmetries ($A_{FB}^{J/\psi}(\theta_{V}),A_{FB}^{J/\psi}(\theta_{V},\theta_{\ell})$) discussed above, where $\chi$ was integrated over the whole range. These observables do not provide any additional information compared to $A_{FB}(\theta_{V})$ and $A_{FB}(\theta_{V},\theta_{\ell})$ discussed before. The 2D scenario in case of $A_{FB}^{J/\psi}(\chi,\theta_{V})$ with Re[$S_{L}=4\,T_{L}$], Im[$S_{L}=4\,T_{L}$] results in about 10-20% deviation from the SM value at low values of $q^{2}$.

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 $\chi$, ($b$) asymmetry depending on $\chi$ and $\theta_{V}$, and ($c$) asymmetry depending on $\chi$, $\theta_{V}$ and $\theta_{\ell}$. The relevant expressions are:

[TABLE]

where $(T_{L}V_{L})^{\ast}={\rm Im}T_{L}{\rm Re}V_{L}-{\rm Im}V_{L}{\rm Re}T_{L},~{}(P\,V_{L})^{\ast}={\rm Im}P\,{\rm Re}V_{L}-{\rm Im}V_{L}{\rm Re}P,~{}(P\,T_{L})^{\ast}={\rm Im}P\,{\rm Re}T_{L}-{\rm Im}T_{L}{\rm Re}P$,

[TABLE]

[TABLE]

The asymmetry $A^{img}_{FB}(\chi,\theta_{V},\theta_{\ell})$ 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 $A^{img}_{FB}(\chi)$ [upper-panel] and $A^{img}_{FB}(\chi,\theta_{V})$ [lower-panel] as a function of $q^{2}$. These observables are only shown for ${\rm Im}S_{L}=4\,{\rm Im}T_{L}$ in the 1D scenario and ${\rm Re}S_{L}=4\,{\rm Re}T_{L}$, ${\rm Im}S_{L}=4\,{\rm Im}T_{L}$ in the 2D scenarios as these were the only cases considered in the global fit in Ref. [43]. The forward-backward asymmetry depending only on $\chi$ 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$ GeV2. The asymmetry $A^{img}_{FB}(\chi,\theta_{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 different NP operators is summarised in Table 6.

[TABLE]

5 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}\to J/\psi l\nu_{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}\rightarrow\eta_{c},J/\psi$ 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 sum rule 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/\psi$ and $\eta_{c}$ are calculated and compared with the results from other approaches. Our predictions for the semileptonic ratios $R_{J/\psi}|_{\rm SM}=0.23\pm 0.01$ and $R_{\eta_{c}}|_{\rm SM}=0.32\pm 0.02$ are in agreement with other derivations and support the existing tension at 2$\sigma$ level with the experiment on $R_{J/\psi}$, Eq. (1). With more data on $B_{c}$ decays from HL-LHC all observables in $B_{c}\to\eta_{c},J/\psi$ 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 $\eta_{c}$ and $J/\psi$ 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(\ast)}$, the $\tau$ and $D^{\ast}$ longitudinal polarization from $B\to D^{\ast}$ decay and the leptonic $B_{c}\to\tau\mu$ branching ratio. We take into account the latest constraints from Ref. [43] and analyse the effects of the NP operators on various observables. The ratio $R_{J/\psi}$ is sensitive to $V_{L}$ in the high $q^{2}$ range whereas $R_{\eta_{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 $\eta_{c}$ decay mode are sensitive to the NP coupling $S_{R}$. The transverse polarization of $\tau$ 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/\psi$, apart from $R_{J/\psi}$.

[TABLE]

In addition, the full 4-fold differential distribution of the decay rate $B_{c}\rightarrow J/\psi\ell\nu_{\ell}$, with $J/\psi$ decaying to a pair of leptons of opposite helicity is considered for the first time in the presence of new physics operators. We find that the asymmetry in the angle $\theta_{V}$ $(A_{FB}^{J/\psi}(\theta_{V}))$ is mostly sensitive to the NP couplings Re,Im[$S_{L}=4T_{L}$], in the 2D NP scenarios. The asymmetries in the angle $\chi$, 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\to D,D^{\ast}$ decays, the asymmetries constructed with $\theta_{V}$, $\chi$ and $(\theta_{V},\chi)$ 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\to(D,D^{\ast})\ell\nu_{\ell}$ decays  [43] can also simultaneously explain the current 2$\sigma$ tension with the experimental $R_{J/\psi}$ ratio. With the extended experimental LHCb program, future studies with more data will be needed to boost or disapprove this evidence of LFU violation in $B_{c}$ decays.

Acknowledgments

This work is partly supported by the EU grant RBI-T-WINNING (grant EU H2020 CSA-2015 number 692194) and by the European Union through the European Regional Development Fund - the Competitiveness and Cohesion Operational Programme (KK.01.1.1.06). M.P. acknowledge support of the Slovenian Research Agency through research core funding No. P1-0035.

Appendix A Form factors calculated in the three-point QCDSR model [4]

We have previously calculated in [4] the same form factors using a more traditional, albeit somewhat modified approach of three-point QCD sum rules (3ptQCDSR) and we present the corresponding results in Table 1. Here we just briefly discuss the method of our calculation and the main difference to the LCSR-inspired approach used in the paper. In 3ptQCDSR mesonic states are interpolated by the currents as

[TABLE]

taken at large virtualities. By inserting a set of hadronic states in the correlation function defined as

[TABLE]

one can extract the form factors by calculating the perturbative part of the correlator and the nonperturbative contributions given in terms of universal vacuum condensates built from the quark and gluon operators of increasing dimension (here we have calculated only the leading nonperturbative contribution coming from the gluon condensate) and matching the QCD result via dispersion relation to a sum over hadronic states. At the end, the expressions are Borel transformed in order to improve the convergence.

Since it is known that in the sum rule calculation of heavy meson decay constants higher orders of perturbation series can contribute as much as 40-50%, whereas the 3-point function is calculated at LO, in order to reduce the uncertainties we have performed the following procedure: in the form factors calculation we have taken for the $s_{0}$ threshold parameters the same values as those that reproduce the corresponding charmonia decay constants obtained from lattice QCD when the decay constants are calculated in the sum rules by taking into account only the LO perturbative part and the gluon-condensate contribution, i.e. with the same approximations as for the form factors, whereas the Borel mass parameter is taken in the region where stability is achieved (we aim at the $\sim 5\%$ stability in the Borel masses in the given Borel window). Furthermore in order to reduce the uncertainties even more, we do not vary the decay constants and thresholds independently, but rather in the 3ptQCDSR calculation we always use the decay constants (varied inside the range allowed by lattice) together with the corresponding thresholds fixed by the decay constants calculation. The hope is that all the higher order/higher dimension operator contributions are then simulated through the appropriate threshold modification in the 3-point QCDSR calculation. The parameters obtained that way are given in Table 8 below. One can notice that in contrast to the LCSR-inspired calculation used in the main text, here the pole mass of the $b$ quark is used, together with the $c$-quark mass derived from the ratio of masses extracted from the lattice calculations.

[TABLE]