Implications of new physics in the decays $B_c \to (J/\psi,\eta_c)\tau\nu$

We study the semileptonic decays of the $B_c$ meson into final charmonium states within the standard model and beyond. The relevant hadronic transition form factors are calculated in the framework of the covariant confined quark model developed by us. We focus on the tau mode of these decays, which may provide some hints of new physics effects. We extend the standard model by assuming a general effective Hamiltonian describing the $b\to c\tau\nu$ transition, which consists of the full set of the four-fermion operators. We then obtain experimental constraints on the Wilson coefficients corresponding to each operator and provide predictions for the branching fractions and other polarization observables in different new physics scenarios.


I. INTRODUCTION
The B c meson is the lowest bound state of two heavy quarks of different flavors, lying below the BD threshold. As a result, while the corresponding cc and bb quarkonia decay strongly and electromagnetically, the B c meson decays weakly, making it possible to study weak decays of doubly heavy mesons. The weak decays of the B c meson proceed via the c-quark decays (∼ 70%), the b-quark decays (∼ 20%), and the weak annihilation (∼ 10%).
Due to its outstanding features, the B c meson and its decays have been studied extensively (for a review, see e.g. Ref. [1] and references therein).
Among many weak decays of the B c meson, the semileptonic decay B c → J/ψℓν has an important meaning. In fact, the first observation of the B c meson by the CDF Collaboration was made in an analysis of this decay [2]. Recently, the LHCb Collaboration has reported their measurement [3] of the ratio of branching fractions mention that very recently, the HPQCD Collaboration has provided their preliminary results for the form factors of the B c → J/ψ and B c → η c transitions using lattice QCD [30].
It should be noted that in our covariant confined quark model (CCQM), all form factors are calculated in the full kinematic range of momentum transfer squared q 2 , making the predictions for physical observables more accurate. In the pQCD approach and QCD sum rules, for example, the form factors are evaluated only for small values of q 2 (large recoil), and then extrapolated to the large q 2 region (small recoil), in which they become less reliable.
In general, the knowledge of the B c → J/ψ(η c ) form factors is much less than that of theB 0 → D ( * ) ones. This is due to, first, the lack of experimental data for the decays B c → J/ψ(η c )ℓν, and second, the appearance of two heavy quark flavors in the initial (bc) and final (cc) states. The latter breaks the heavy quark symmetry (HQS), leaving the residual heavy quark spin symmetry (HQSS), which allows reducing the number of form factors in the infinite heavy quark limit [19,31,32]. However, the HQSS does not fix the normalization of the form factors as the HQS does, for example, in the case ofB 0 → D ( * ) .
The possible NP effects in the semileptonic decays B c → J/ψ(η c )τ ν have been discussed recently in several papers [33][34][35][36][37][38][39][40]. As what has been done in the studies of the R D ( * ) anomalies, one can choose between a specific-model approach, such as charged Higgs models, leptoquark models etc., or a model-independent approach based on a general effective Hamiltonian describing the b → cτ ν transition. In this paper we adopt the second approach by proposing an SM-extended effective Hamiltonian consisting of the full set of the fourfermion operators. Constraints on the corresponding Wilson coefficients are obtained from the experimental data for the ratios R J/ψ and R D ( * ) , as well as the LEP1 data, which requires B(B c → τ ν) ≤ 10% [41]. Another useful constraint is provided by using the lifetime of B c as discussed in Ref. [42]. However, in this paper we only use the constraint from the LEP1 data, which is more stringent than the latter. We then analyze the effects of these NP operators on several physical observables including the ratios of branching fractions R J/ψ(ηc) , the forward-backward asymmetries, the convexity parameter, and the polarization components of the τ in the final state. We also provide our predictions for these physical observables in the SM and in the presence of NP.
The paper is organized as follows. In Sec. II we introduce the general b → cℓν effective Hamiltonian and parametrize the hadronic matrix elements in terms of the invariant form factors. We then obtain the decay distributions in the presence of NP operators using the helicity technique. In Sec. III we present our result for the form factors in the whole q 2 range. A detailed comparison of the form factors calculated in the CCQM with those in other approaches is also provided. In Sec. IV we obtain constraints on the NP Wilson coefficients from available experimental data. Theoretical predictions for the physical observables in the SM and beyond are presented in Sec. V. Finally, we briefly conclude in Sec. VI.

II. EFFECTIVE HAMILTONIAN, HELICITY AMPLITUDES, AND DECAY DIS-TRIBUTION
In the model-independent approach, NP effects are introduced explicitly by proposing an effective Hamiltonian for the weak decays that includes both SM and beyond-SM contributions. In this study, the general effective Hamiltonian for the quark-level transition b → cℓν (ℓ = e, µ, τ ) is given by (i = L, R) where the four-fermion operators O X are defined as Here, σ µν = i [γ µ , γ ν ] /2, P L,R = (1 ∓ γ 5 )/2 are the left and right projection operators, and X's are the complex Wilson coefficients characterizing the NP contributions. The tensor operator with right-handed quark current is identically equal to zero and is therefore omitted.
In the SM one has V L,R = S L,R = T L = 0. We have assumed that neutrinos are left-handed.
Besides, the delta function in Eq. (2) implies that NP effects are supposed to appear in the tau mode only. The proposed Hamiltonian can be considered as a natural way to go beyond the SM since it is generalized from the well established SM Hamiltonian with the V − A structure by adding more currents. One may also consider right-handed neutrinos and may as well assume that NP appears in all lepton generations. However, current experimental data suggest that NP effects in the case of light leptons (if any) are very small. A recent discussion of these NP operators and their possible appearance in the light lepton modes can be found in Ref. [43].
Starting with the effective Hamiltonian, one writes down the matrix element of the semileptonic decays B c → J/ψ(η c )τ ν, which has the form where Γ X is the Dirac matrix corresponding to the operator O X . The hadronic part in the matrix element is parametrized by a set of invariant form factors depending on the momentum transfer squared q 2 between the two hadrons as follows: where P = p 1 + p 2 , q = p 1 − p 2 , and ǫ 2 is the polarization vector of the J/ψ meson which satisfies the condition ǫ † 2 · p 2 = 0. The particles are on their mass shells: p 2 1 = m 2 1 = m 2 Bc and p 2 2 = m 2 2 = m 2 J/ψ(ηc) . We define a polar angle θ as the angle between q = p 1 − p 2 and the three-momentum of the charged lepton in the (ℓν ℓ ) rest frame. The angular decay distribution then reads where |p 2 | = λ 1/2 (m 2 1 , m 2 2 , q 2 )/2m 1 is the momentum of the daughter meson in the B c rest frame, and H µν L µν is the contraction of hadron and lepton tensors. The covariant contraction H µν L µν can be converted to a sum of bilinear products of hadronic and leptonic helicity amplitudes using the completeness relation for the polarization four-vectors of the process [44]. This technique is known as the helicity technique, which has been described in great detail in our previous papers [44][45][46][47]. In Ref. [47] we have shown how to acquire the decay distribution for the semileptonic decaysB 0 → D ( * ) τ ν in the presence of NP operators and provided a full description of the helicity amplitudes, which can be applied to the case of the decays B c → J/ψ(η c )τ ν. Therefore, we find no reason to repeat the procedure in this paper. However, for completeness, we present here the final result for the decay distributions.
The angular distribution for the decay B c → η c τ ν is written as follows: where  7). Their explicit expressions are presented in Ref. [47]. Note that we do not consider interference terms between different NP operators since we assume the dominance of only one NP operator besides the SM contribution. The corresponding distribution for the decay B c → J/ψτ ν is rather cumbersome and therefore is not shown here. One can find it in Appendix C of Ref. [47].
After integrating the angular distribution over cos θ one has In this paper, we also impose the constraint from the leptonic decay channel of B c on the Wilson coefficients. Therefore we present here the leptonic branching in the presence of NP operators. In the SM, the purely leptonic decays B c → ℓν proceed via the annihilation of the quark pair into an off shell W boson. Assuming the effective Hamiltonian Eq. (2), the tau mode of these decays receives NP contributions from all operators except O T L . The branching fraction of the leptonic decay in the presence of NP is given by [48] where of B c , and f P Bc is a new constant corresponding to the new quark current structure. One has In the CCQM, we obtain the following values for these constants (all in MeV):

III. FORM FACTORS IN THE COVARIANT CONFINED QUARK MODEL
The CCQM is an effective quantum field approach to hadron physics, which is based on a relativistic invariant Lagrangian describing the interaction of a hadron with its constituent quarks (see e.g. Refs. [49][50][51][52][53][54][55]). The hadron is described by a field H(x), which satisfies the corresponding equation of motion, while the quark part is introduced by an interpolating quark current J H (x) with the hadron quantum numbers. In the case of mesons, the Lagrangian is written as where g H is the quark-meson coupling, Γ H is the Dirac matrix ensuring the quantum numbers of the meson, and the so-called vertex function F H effectively describes the quark distribution inside the meson. From the requirement for the translational invariance of F H , we adopt the The coupling g H is determined by using the so-called compositeness condition [56], which imposes that the wave function renormalization constant of the hadron is equal to zero Z H = 0. For mesons, the condition has the form is the derivative of the hadron mass operator, which corresponds to the self-energy diagram in Fig. 1 and has the following form for pseudoscalar and vector mesons, respectively. Here, S 1,2 are quark propagators, for which we use the Fock-Schwinger representation It should be noted that all loop integrations are carried out in Euclidean space.
Similarly to the hadron mass operator, matrix elements of hadronic transitions are represented by quark-loop diagrams, which are described as convolutions of the corresponding quark propagators and vertex functions. Using various techniques described in our previous papers, any hadronic matrix element Π can be finally written in the form where F is the resulting integrand corresponding to a given diagram. It is more convenient to turn the set of Fock-Schwinger parameters into a simplex by α i as follows: The integral in Eq. (19) begins to diverge when t → ∞, if the kinematic variables allow the appearance of branching point corresponding to the creation of free quarks. However, these possible threshold singularities disappear if one cuts off the integral at the upper limit: The parameter λ effectively guarantees the confinement of quarks inside a hadron and is called the infrared cutoff parameter.
The parameters of the B c → J/ψ(η c ) form factors are listed in Table I. Their q 2 dependence in the full momentum transfer range 0 ≤ q 2 ≤ q 2 max = (m Bc − m J/ψ(ηc) ) 2 is shown in Fig. 2. Firstly, we focus on those form factors that are needed to describe the B c → J/ψ(η c ) transitions within the SM (without any NP operators), namely, F ± , A 0,± , and V . It is worth noting that all these form factors have a pronounced (q 2 ) −2 contribution (the ratio b/a lies between 0.14 and 0.50) in comparison with the caseB 0 → D ( * ) , where all form factors (except for A 0 ) have a very small ratio b/a ∼ 0.05−0.08, and therefore show a monopolelike behavior [46].  For easy comparison, we relate all form factors to the well-known Bauer-Stech-Wirbel form factors [57], namely, F +,0 for B c → η c , and A 0,1,2 and V for B c → J/ψ. Note that in Ref. [57] the notation F 1 was used instead of F + . In Fig. 3 and Fig. 4 we compare our form factors with those obtained in other approaches, namely, perturbative QCD [25], QCD sum rules (QCDSR) [21], the Ebert-Faustov-Galkin relativistic quark model [22], the Hernandez-Nieves-Verde-Velasco (HNV) nonrelativistic quark model [23], and the covariant light-front quark model (CLFQM) [18]. It is interesting to note that our form factors are very close to those computed in the CLFQM [18].  Using the heavy quark spin symmetry, the authors of Ref. [19] have obtained several relations between the form factors of the B c → J/ψ(η c ) transitions. In particular, the relation between the form factors F + and F − can be used to prove the linear behavior of the ratio F 0 (q 2 )/F + (q 2 ), where the slope α only depends on the masses of the involved quarks and hadrons. We find α = 0.020 GeV −2 from the numerical values in Ref. [19]. Similarly, we obtain α = 0.018 GeV −2 from results of Refs. [21,25], α = 0.021 GeV −2 from Ref. [18]. However, Ref. [22] and Ref. [23] yield much smaller values, which are α = 0.005 GeV −2 and α = 0.007 GeV −2 , respectively. In our model, the ratio F 0 (q 2 )/F + (q 2 ) exhibits an almost linear behavior in the whole q 2 range as demonstrated in Fig. 5, from which we obtain α = 0.017 GeV −2 . The value of the slope α plays an important role in studying the shape of the form factors, which can be determined more accurately by future lattice calculations. It is also worth mentioning the very recent lattice results for the B c → J/ψ form factors provided by the HPQCD Collaboration [30]. In this study, they found A 1 (0) = 0.49, A 1 (q 2 max ) = 0.79, and V (0) = 0.77, which are very close to our values A 1 (0) = 0.56, A 1 (q 2 max ) = 0.79, and V (0) = 0.78. Almost all the recent studies on possible NP in the decays B c → (J/ψ, η c )τ ν employ the form factors F 0,+ , A 0,1,2 , and V calculated in pQCD approach [25]. The remaining form factors corresponding to the NP operators are obtained by using the quark-level equations of motion (EOMs). In this paper we provide the full set of form factors in the SM as well as in the presence of NP operators without relying on the EOMs. However, this does not mean that our form factors do not satisfy the EOMs. A brief discussion of the EOMs in our model can be found in Ref. [47]. Our form factors therefore can be used to analyze NP effects in the decays B c → (J/ψ, η c )τ ν in a self-consistent manner, and independently from other studies.

IV. EXPERIMENTAL CONSTRAINTS
Constraints on the Wilson coefficients appearing in the effective Hamiltonian Eq. (2) are obtained by using experimental data for the ratios of branching fractions R D = 0.407±0.046, R D * = 0.304 ± 0.015 [13], and R J/ψ = 0.71 ± 0.25 [3], as well as the requirement B(B c → τ ν) ≤ 10% from the LEP1 data [41]. It should be mentioned that within the SM our calculation yields R D = 0.267, R D * = 0.238, and R J/ψ = 0.24. We take into account a theoretical error of 10% for our ratios. Besides, we assume the dominance of only one NP operator besides the SM contribution, which means that only one NP Wilson coefficient is considered at a time. In Fig. 6 we show the constraints on the scalar Wilson coefficients S L,R within 2σ. It is seen that the recent experimental value of R J/ψ does not give any additional constraint on S L,R to what have been obtained by using R D ( * ) . In particular, S R is excluded within 2σ using only R D ( * ) . However, in the case of S L , the constraint from B(B c → τ ν) plays the main role in ruling out S L . In general, the branching of the tauonic B c decay imposes a severe constraint on the scalar NP scenarios. Many models of NP involving new particles, such as charged Higgses or leptoquarks, also suffer from the same constraint, and therefore need additional modifications to accommodate the current experimental data (see e.g. Refs. [58][59][60]). In the upper panels of Fig. 7 we present the constraints on the vector V L,R and tensor T L Wilson coefficients. There is no available space for these coefficients within 1σ. Moreover, they are excluded mainly due to the additional constraint from R J/ψ , rather than from B(B c → τ ν). This holds exactly in the case of T L since the operator O T L has no effect on B(B c → τ ν). In the lower panels of Fig. 7 we show the allowed regions for V L,R and T L within 2σ. In each allowed region at 2σ we find a best-fit value for each NP coupling. The best-fit couplings read V L = −1.05 + i1.15, V R = 0.04 + i0.60, T L = 0.38 − i0.06, and are marked with an asterisk.

V. THEORETICAL PREDICTIONS
In this section we use the 2σ allowed regions for V L,R and T L obtained in the previous section to analyze their effects on several physical observables. Firstly, in Fig. 8 we show the q 2 dependence of the ratios R J/ψ and R ηc in different NP scenarios. It is obvious that all the NP operators O V L , O V R , and O T L increase the ratios. However, it is interesting to note that O T L can change the shape of R J/ψ (q 2 ) and may imply a peak in the distribution. This unique behavior can help identify the tensor origin of NP by studying the q 2 distribution of the decay B c → J/ψτ ν. The average values of the ratios R J/ψ and R ηc over the whole q 2 region are given in Table II. The row labeled by SM contains our predictions within the SM using our form factors. The predicted ranges for the ratios in the presence of NP are given in correspondence with the 2σ allowed regions of the NP couplings shown in Fig. 7. Here, the most visible effect comes from the operator O V R , which can increase the average ratio < R ηc > by a factor of 2.
Next, we consider the polarization observables in these decays. For this purpose we write the differential (q 2 , cos θ) distribution as where W (θ) is the polar angular distribution, which is described by a tilted parabola. For convenience we define a normalized polar angular distribution W (θ) as follows: The normalized angular decay distribution W (θ) obviously integrates to 1 after cos θ integration. The linear coefficient b/2(a+c/3) can be projected out by defining a forward-backward asymmetry given by The quadratic coefficient c/2(a + c/3) is obtained by taking the second derivative of W (θ).
We therefore define a convexity parameter by writing In the upper panels of Fig. 9 we present the q 2 dependence of the forward-backward In the lower panels of Fig. 9 we show the convexity parameter C τ F (q 2 ). It is seen that the operator O V R has a very small effect on C τ F , and only in the case of B c → J/ψ. In contrast to this, C τ F is extremely sensitive to the tensor operator O T L . In particular, O T L can change where s µ i are the polarization four-vectors of the τ in the W − rest frame. One has Here, p τ and p 2 are the three-momenta of the τ and the final meson (J/ψ or η c ), respectively, in the W − rest frame. A detailed analysis of the tau polarization with the help of its subsequent decays can be found in Refs. [63][64][65].
The q 2 dependence of the tau polarizations are presented in Fig. 10.  Table III with the same notations as for Table II.   namely, the ratios R J/ψ (q 2 ) and R ηc (q 2 ), the forward-backward asymmetry A F B (q 2 ), the convexity parameter C τ F (q 2 ), and the polarizations of the τ in the final state. Some of the effects may help distinguish between NP operators. We have also provided predictions for the q 2 average of the mentioned observables, which will be useful for other theoretical studies and future experiments.