The Study of Rare $B_c\rightarrow D^{(*)}_{s,d}l\bar{l}$ Decays

In this paper, we study rare decays $B_c\rightarrow D^{(*)}_{s,d}l\bar{l}$ within the Standard Model. The penguin, box, annihilation, color-favored cascade and color-suppressed cascade contributions are included. Based on our calculation, the annihilation and color-favored cascade diagrams play important roles in the differential branching fractions, forward-backward asymmetries, longitudinal polarizations of the final vector mesons and leptonic longitudinal polarization asymmetries. More importantly, color-favored cascade decays largely enhance the resonance cascade contributions. To avoid the resonance cascade contribution pollution, new cutting regions are put forward.

mesons is expected to be 5 × 10 10 [2,3] per year when it runs at the highest luminosity 10 34 cm −2 s −1 and highest energy 14 TeV. So most B c decay channels will be observed with this huge number. It is believed that the predictions about the B c meson in theory will be checked by the experiment data in the near future.
Among the channels of the B c meson transitions, the rare decays b → s(d)ll and b → s(d)γ have been emphasized recently. These decays are the single-quark flavor-changing neutral current (FCNC) processes. Based on the Glashow-Iliopoulos-Maiani (GIM) [4] mechanism in the Standard Model (SM), they are forbidden at tree level but induced by electro-weak loop diagrams. It implies that the SM contribution is greatly suppressed. And, as the SM contribution is suppressed, the new physics (NP) effects may become important. So studying the rare decays will be helpful to test the SM stringently and detect the NP effects indirectly. Among these rare decay channels, the semi-leptonic processes b → s(d)ll are favored because of their abundant observables, for instance, the differential branching fraction, the forward-backward asymmetry, the longitudinal polarization of the final vector meson, the leptonic longitudinal polarization asymmetry and so on. More information on the Wilson coefficients of the rare decays can be extracted from these observables.
In existing studies, the b → sll processes have been discussed widely through the B → K ( * ) ll decays. In the B → K ( * ) ll processes, the b → sll contribution is dominant, while the annihilation diagrams are CKM suppressed |V * ub V us |/|V * ts V tb | ∼ λ 2 and the spectator scattering effect is at the next α s order. As shown in [5][6][7], by including only the b → sll contribution in the SM, s,d ll rare decays are CKM allowed |V * cb V cs(d) |/|V * ts(d) V tb | ∼ 1 and enhanced by a 3 times larger color factor. Thus, the annihilation contributions should be considered seriously. In this way, the study of the B c → D ( * ) s,d ll rare decays provides outstanding features. First, the NP effects will be more distinct by considering the annihilation processes. As pointed out in [8,9], with the help of the interference between the b → s(d)ll contribution and the SM annihilation contribution, the predictions of SM with a fourth generation (SM4) and the Super-Symmetric (SUSY) Models will deviate stronger from the SM predictions in the processes B c → D * s ll than in the processes B → K * ll. Furthermore, in B c rare decays, the particles beyond SM may contribute not only to the b → sll transitions but also to the annihilation diagrams. If so, this double-contribution mechanism may make the NP effect more evident. Second, the processes B c → D s,d ll can help us comprehend the large isospin asymmetry phenomenon in the B → Kµμ decay. Recently, LHCb [10] has reported that the isospin asymmetry A I of the B → K * µμ process is in agreement with the naive SM expectation of a vanishing asymmetry. However, in the B → Kµμ process, there is a 4σ deviation from zero for A I . Even if calculations [11][12][13] including the annihilation and spectator-scattering diagrams are carried out, the isospin asymmetry of the process B → Kµμ is also less than 0.05. One may infer that studying the spectator effects should be emphasized.
In the decays B c → D s,d ll, spectator effects will be obvious based on the particular annihilation mechanism. So the rare B c decays can offer a new and helpful field to reveal the reason of the large isospin asymmetry in the B → Kµμ decay.
In previous works [5,[14][15][16], which include the b → s(d)ll short-distance contribution and the b →ccs(d) → s(d)ll long-distance diagrams, the B c → D s,d ll and B c → D * s,d ll processes have been calculated, while in Refs. [17][18][19] only the b → s(d)ll short-distance contribution is considered. Recently, using the annihilation form-factors of B c → D * s γ [20][21][22], the B c → D * s ll transitions combined with the annihilation effects have been analysed in Refs. [8,9,23]. Ac- s,d ll annihilation diagrams are considered, more form-factors will appear. In the aspect of the long-distance interaction, the papers [5,8,[14][15][16] have calculated the colorsuppressed cascade diagrams. The color-suppressed cascade long-distance contribution dominates the B → K ( * ) J/ψ(ψ(2S)) → K ( * )l l processes. However, both color-suppressed and color-favored long-distance diagrams contribute to B c → D To investigate the B c → D ( * ) s,d ll decays, the Operator Production Expansion (OPE) mechanism and factorization ansatz can be employed. In this method, the amplitude can be separated into two parts: the short-distance Wilson coefficients (If the decays involve the resonance cascade processes, extra long-distance terms will appear.) and the long-distance hadronic matrix elements which are operators sandwiched by the initial and the final states. In the framework of Re-normalization Group Equations (RGE), the Wilson coefficients can be obtained at next-to-leading (NL) order or next-to-next-to-leading (NNL) order QCD corrections [25,26] perturbatively. But calculating the hadronic matrix elements is a non-perturbative problem and a model dependent method has to be chosen to do it.
Till now, the hadronic matrix elements have been investigated in several approaches: the relativistic constituent quark model [14,15], the light-cone quark model [15][16][17], the three point QCD sum rules [18,19] and the QCD-motivated relativistic quark model [5]. In this paper, we choose the Mandelstam Formalism (MF) [27] and the Bethe-Salpeter (BS) equation [28,29] to investigate the hadronic matrix elements. This method has several particular features. First, the BS equation, based on the quantum field theory, is a relativistic equation to describe a two-body bound state. Using this method, the Gaussian wave function is abandoned and a relativistic form of the wave function has been introduced [30][31][32]. In this way, the wave function has the same parity and charge parity as the bound state. Different forms of wave functions denote different states, e.g., the wave function of a pseudoscalar meson is much different from the one of a vector meson. Only in the non-relativistic limit, they are similar [33]. Second, in the decays s,d ll, the mass of initial meson B c is much larger than that of the final meson D ( * ) s,d . As a result, the final meson recoil can be large and the relativistic effect should be taken seriously.
Based on the Mandelstam Formalism, our method keeps the relativistic effect not only from the relativistic wave functions but also from the kinematics [34]. At last, the weak annihilation amplitude involving the heavy-light meson B can be calculated within the Heavy Quark Effective Theory (HQET) [35,36] under the m u(d) /m b ∼ 0 limit. However, unlike the B and D meson, the B c meson consists of two heavy quarks and the relationship m c /m b ∼ m u(d) /m b ∼ 0 is not suitable. So in this paper, we calculate the annihilation hadronic matrix elements based on the BS method. In our calculation, dynamics of both heavy quarks is considered. And our method can be used to deal with not only the double-heavy quark meson but also the heavylight meson. Besides, the weak-annihilation hadronic currents satisfy the gauge-invariance at the (P i − P f ) 2 = 0 GeV 2 . This paper is organized as follows. In section 2, the theoretical details of the effective hamiltonian, hadronic transition matrix elements and the observables of the B c → D ( * ) s,d ll processes are given. In section 3, the numerical results and discussions are presented. In section 4, we summarize this paper. We have the effective Hamiltonian for annihilation, color-suppressed cascade and color-favored cascade diagrams: where G F is the Fermi constant and C 1 , C 2 are willson coefficients. In this paper, when the annihilation, color-suppressed and color-favored cascade diagrams are calculated, the effects of Q 3−6 are neglected because of their small Wilson coefficients.
The penguin and box diagrams' effective Hamiltonian, that is, the b → s(d)ll effective Hamiltonian, is given by: where V cb , V cs(d) , V tb and V ts(d) are the CKM matrix elements. The set of local operators is [14,25]: For b → s(d)ll processes, with the help of RGE [37] approach, C 7,9 (m b ) are linear combinations of according C 7,9 (M W ) and C 1−6 (M W ). The operators Q 1−6 themselves have matrix elements that contribute to b → sll and the perturbative parts can be included into C ef f 7,9 . In our calculation, considering that the effect of Q 8 belongs to the higher α s order, we do not include its contribution.
Thus, the amplitude induced by effective Hamiltonian is shown as: In the above equation, we have defined the hadronic matrix elements W µ and W T µ : where P i and P f are momenta of the initial and final mesons, respectively. The antisymmetric tensor is defined as σ µν = i 2 [γ µ , γ ν ] and m b is the mass of b quark.
The effective Wilson coefficients C ef f 7 and C ef f 9 are defined as: where h(m c , s) (m c = m c /m b ), h(1, s) and h(0, s) describe the contributions of the four-fermi operator Q 1−6 loops. In Refs [14,37,38], with the help of RGE, the Wilson coefficients C i (i = 1 − 10) are evolved from M W scale to the m b scale. In this paper, we adopt the same expressions of h functions and the same numerical values of the Wilson coefficients as those in Ref. [14].
For the annihilation diagrams, according to the naive factorization, we write the amplitude as: where in the second equality the Fierz-arrangement identity is employed. The annihilation hadronic matrix elements can be expressed as: In the above equation, W µ 1,2,3,4 are defined as: where p q 1(2−4) and m q 1(2−4) are momenta and masses of the propagated quarks, respectively.
The LD processes considered in this paper are those that are induced by the resonance s(d)l l, whose contributions can be described by the relationship The relationship written here is not used to compute any observables in this paper but just employed to make the following sentences transparent.) The resonances V denote J P C = 1 −− mesons which could be theūu,dd,ss andcc bound states. During our calculation, we ignore the effects s(d)l l cascade decays. On one hand, the strong decays allowed by Okubo-Zweig-Iizuka (OZI) rules are the dominant channels of the ρ, ω and φ mesons, while the strong decays of J/ψ(ψ(2S)) are suppressed by OZI rules. So the processes ρ(ω, φ) → l + l − which are induced by electromagnetic interaction have much smaller branching fractions than J/ψ(ψ(2S)) → l + l − . On the other hand, considering the small CKM matrix elements V ub and V us , the branching fractions of B c → D ( * ) s ρ(ω) processes will be suppressed. And because of s(d)l l processes in this paper.
As pointed out in Refs. [39,40], the propagator of V can be written in the Breit-Wigner form. So, from the naive factorization, the amplitude can be written as: where J are the mass and width of the resonance meson, respectively. From Fierz-arrangement identity and naive factorization, we write the where on the right side of the above equation, the first term is from the CS diagram and the second term gives the CF contribution. If only the CS and b → s(d)ll contributions are included, our amplitude M b→s(d)l + l − + M LD will be similar with those in Refs. [5,14,[41][42][43]. Taking account of the unitarity condition on the elastic Breit-Wigner resonance amplitude [42], we take ϕ = 0 in this paper.
To simplify the calculation, we define: where in the second step the relationship Γ(V →ll) = πα 2 em 16f 2 V /(27M V ) is used. And the C CS 9 is given by: Similarly, we also define: By including the BP, Ann, CF cascade and CS cascade diagrams, the amplitude calculated in this paper can be written as: Considering the Lorentz and parity symmetries, W µ and W T µ could be written in terms of the form factors: The form-factors of annihilation and color-favored diagrams are defined as: where and V CF are the form factors.

Hadronic Transition Matrix Elements in the Bethe-Salpeter Method
In this subsection, we show the details of how to calculate the hadronic matrix elements in Eqs. (12,13). Within the Bethe-Salpeter method, a meson is considered as a bound state of the quark and the anti-quark. For a meson with mass M , momentum P and relative momentum q, we can define: where p 1 (p 2 ) and m 1 (m 2 ) are the momentum and mass of the constituent quark(anti-quark), respectively. To describe the bound state, the wave function is essential. In this paper, we obtain it by solving the BS equation: where χ(q) is the BS wave function and V (P, k, q) is the interaction kernel.
To solve the equation, we need to reduce it to its instantaneous version, the Salpeter equation [30]. Then we need a proper interaction kernel and the forms of the wave functions. In our method, the Cornell potential, which is one of the most successful and widely used potentials in describing heavy mesons, is adopted as the interaction kernel.
The behavior of a bound state is determined by its J P quantum number, so we give the form of the wave function based on the meson's J P quantum number. As a result, the wave functions of the different mesons fulfill different BS equations [30,33]. For a pseudoscalar meson, J P = 0 − , the positive energy wave function can be written as [30]: One can check that the wave function has the correct J P quantum number. In Eq. (14), the parameters are defined as: where ω i is defined as ) and ϕ 2 (q 2 P ⊥ ) are obtained by solving the full Salpeter equations [30].
For vector mesons, J P = 1 − , the positive energy wave function can be written as [31]: where the coefficients are shown as: In Eqs. (16,17), the same definitions of P , M , ω i , m i and q P ⊥ as for the pseudo-scalar are used.
For J/ψ and ψ(2S) mesons, i = 1 denotes the c quark and i = 2 denotes thec quark. For D * − s(d) , i = 1 denotes s(d) quark and i = 2 denotesc quark. The numerical values of wave functions ϕ 3 (q P ⊥ ), ϕ 4 (q P ⊥ ), ϕ 5 (q P ⊥ ) and ϕ 6 (q P ⊥ ) can be obtained by solving the full Salpeter equations for a vector meson. And the details can be found in Ref. [31].
Using the MF [27], the hadronic matrix elements W µ and W µ T can be expressed as an overlapping integral over the initial and final meson wave functions [33]. We find that the contribution from the positive energy wave function is dominant. Thus, the other parts have been ignored. After integrating over q 0 P ⊥ , the hadronic matrix elements can be obtained: where P i and M i are the momentum and mass of the initial meson, respectively. ϕ ++ i and ϕ ++ f are the positive energy wave functions of the initial and final mesons, respectively.
Using this relationship, the calculations can be highly simplified but a tiny deviation will emerge.
For the color-favored long-distance hadronic matrix element W µ CF as shown in Eq. (11), it is proportional to the product of the transition matrix element V |cγ ν (1 − γ 5 )b|B c multiplied by decay constant D

s,d ll processes
As pointed out in Refs. [5,14], during the calculation of the observables, it is convenient to project the hadronic matrix elements to the 4-Helicity-Basis ε †µ H (±, 0, t). With this projection, both the hadronic and leptonic currents can be calculated in their own C.M.S separately. And the leptonic angular distributions can be easily expressed in the terms of helicity amplitudes.
The helicity amplitudes in this paper are shown as: where where because the final meson D s(d) is a pseudo-scalar meson and has no polarization direction, the transverse helicity amplitudes of the B c → D s(d) l ± l ∓ processes are 0 (Of course, if new physics operators contribute, they will be in a different situation.).
Based on the calculations in Ref. [14], the differential branching fraction is: where M H is defined as: Besides that, we also compute the forward-backward asymmetries A F B and the longitudinal polarizations P L of the final vector mesons in the decays B c → D * s,dl l. In the investigation of the B → K * l l processes, A F B and P L have been widely paid attention to, both theoretically and experimentally. We hope to obtain more information on the Wilson coefficients through studying these observables. From the results in Ref. [14], we have: In this paper, only the longitudinal polarizations P L of the final vector mesons are investigated.
The transverse polarizations P T of the final vector mesons can be obtained from the relationship P T = 1 − P L obviously. Furthermore, we also study the leptonic longitudinal polarization asymmetry A LP L , which is defined as: where h = +1(−1) denotes right (left) handed l − . And the second step followed the derivation in Ref [14].

Numerical Results and Discussions
In this section, the parameters, the form factors, the differential branching fractions, forwardbackward asymmetries, longitudinal polarizations of the final vector mesons and leptonic longitudinal polarization asymmetries are presented.  [44], the mass spectra with these input constituent quark masses are in good agreement with the experimental data. Thus, this set of numerical values is adopted in this paper. The numerical values of Cornell-potential-parameters are adopted as
Moreover, the masses and the lifetimes of B c , D ( * ) s(d) , J/ψ, ψ(2S) are taken from Particle Data Group (PDG) [45], as the values of α em , G F and V CKM .
In this paper, the theoretical uncertainties are estimated including two aspects. On one hand, theoretical uncertainties caused by BS-model dependent parameters are studied. This kind of errors is calculated by changing the numerical values of our model dependent parameters by ±5%, which include the masses of constituent quarks and the Cornell-potential-parameters. We find that, in BS method, the observables are more sensitive to Λ BS and λ BS than the other BS-model dependent parameters.
On the other hand, the systematic uncertainties which the naive factorisation hypothesis arouses are investigated. In the investigation of B → K ( * )l l processes, to include non-factorisable contributions of LD cascade processes B → K ( * ) J/ψ(ψ(2S)) → K ( * )l l, the κ factor which is determined by comparing the experimental data with the theoretical results is introduced frequently [5,14,46]. However, the situation in B c → D ( * ) s,d J/ψ(ψ(2S)) → D ( * ) s,dl l is different. The cascade processes B → K ( * ) J/ψ(ψ(2S)) → K ( * )l l include only CS contributions and the κ factor can be used to include the non-factorisable effects of CS diagrams, but transitions s,dl l contain both CS and CF diagrams. Thus, another way is adopted in this paper. In the previous calculations on the non-leptonic decays of B meson [47][48][49][50], to contain the non-factorisable contributions, the number of colors N c in the expression Considering that C ef f 7 has a small value, we can write the relationship H s(d) has a large recoil momentum. As a result, the overlapping region of initial and final wave functions is small and the form factors of W µ and W T µ will be suppressed around Q 2 ∼ 0 GeV 2 . In the high Q 2 region, the final meson D −( * ) s(d) has a tiny recoil momentum and the overlapping region of initial and final wave functions will be enhanced. So we find that the form factors of W µ and W T µ become larger with Q 2 increasing.
However, the calculation of W µ ann is quite different. As shown in Eqs. (7,20), it is a product of the initial (final) meson decay constant multiplied by the final (initial) annihilation hadronic matrix element. The initial and final wave functions are integrated independently. To do the integral, we use the following relationship, In this way, the form-factors of W µ ann are complex (The same scenario has also been reported by Ref. [51] in analysis of the t-channel annihilation diagrams in the B 0 → D 0 µμ process by pQCD method.) and we define them as:

Observables and Discussions
With the form-factors above, the differential branching fractions dBr/dQ 2 , forward-backward asymmetries A F B , leptonic longitudinal polarization asymmetries A LP L and longitudinal polarizations P L of the final vector mesons are calculated using Eqs. (21)(22)(23)(24)(25) and shown in Figs. 6-13. In our calculation, if only Q 7 and Ann contributions are considered, we find that the relationship dBr Q 7 +Ann /dQ 2 > dBr Ann /dQ 2 > dBr Q 7 /dQ 2 is established in the Q 2 ∼ 0 GeV 2 region. In calculation of B c → D * s,d γ [20][21][22], the same relationship also holds.
We 12 (e, f), P L are almost proportional to Q 2 within the low Q 2 region but inversely related to Q 2 in the high Q 2 area. In the Q 2 ∼ Q 2 M AX region, the final meson has a tiny recoil momentum and the polarization is almost averaged. So P L ∼ 1 3 is found near the Q 2 M AX point. In the Q 2 ∼ 0 GeV 2 region, the Q 7 and annihilation diagrams play important parts in P L . Considering that the real γ has no longitudinal polarization, in a physical sense, P L should be towards 0 nearby the Q 2 = 0 GeV 2 point. This feature is in agreement with the one in Ref. [5] but has obvious difference from those in Refs. [8,23]. Beside that, as shown in Fig. 10 (c) and Fig. 12 (c), P L with only PB processes differ slightly from those which include both PB and CS contributions.
In Figs. 6-13 (a, b), the differential branching fractions of B c → D ( * ) s(d)l l are plotted. From the solid lines of Fig. 6 (b), Fig. 8 (b), Fig. 10 (b) and Fig. 12 (b), it is found that the differential branching fractions near the resonance region of J/ψ are a lot larger than the ones around the resonance region of ψ(2S). The cases of the Fig. 6 (a), Fig. 8 (a), Fig. 10 (a) and Fig. 12 (a)  s(d)l l are calculated, the cutting regions on the di-muon mass around the J/ψ and ψ(2S) masses are required. In Refs. [15,16], through analysing the contributions of the PB and CS cascade diagrams, the experimental cutting regions were defined. However, based on our differential branching fractions in Figs. 6-13 (a, b), the CF cascade diagrams affect the differential branching fractions substantially. So new experimental cutting regions will be necessary.
In this paper, we attempt to define the experimental cutting regions: 5 GeV 2 < Q 2 < 15 GeV 2 and the regions of interest are given as: Region (1) : 0.2 GeV 2 < Q 2 < 5 GeV 2 , In

Conclusion
In summery, the processes B c → D