Angular distribution of the FCNC process $B_{c}\to{D}_{s}^{*}(\to{D}_{s}\pi)\ell^{+}\ell^{-}$

In this work, we study the flavor-changing neutral-current process $B_{c}\to{D}_{s}^{*}(\to{D}_{s}\pi)\ell^{+}\ell^{-}$ ($\ell$= $e$, $\mu$, $\tau$). The relevant weak transition form factors are obtained by using the covariant light-front quark model, in which, the main inputs, i.e., the meson wave functions of $B_{c}$ and $D_{s}^{*}$, are adopted as the numerical wave functions from the solution of the Schr\"{o}dinger equation with the modified Godfrey-Isgur model. With the obtained form factors, we further investigate the relevant branching fractions and their ratios, and some angular observables, i.e., the forward-backward asymmetry $A_{FB}$, the polarization fractions $F_{L(T)}$, and the $CP$-averaged angular coefficients $S_{i}$ and the $CP$ asymmetry coefficients $A_{i}$. We also present our results of the clean angular observables $P_{1,2,3}$ and $P^{\prime}_{4,5,6,8}$, which can reduce the uncertainties from the form factors. Our results show that the corresponding branching fractions of the electron or muon channels can reach up to $10^{-8}$. With more data being accumulated in the LHCb experiment, our results are helpful for exploring this process, and deepen our understanding of the physics around the $b\to{s}\ell^{+}\ell^{-}$ process.


I. INTRODUCTION
The flavor-changing neutral-current (FCNC) process, like the b → sℓ + ℓ − (ℓ=e, µ, τ) we are concerned with has attracted the attention of both theorists and experimentalists, and of course has been widely studied.The FCNC process is forbidden at the tree level, and can only operate through loop diagrams in the Standard Model (SM).At the lowest order, three amplitudes contribute to the decay width, i.e., the photo penguin diagram, the Z penguin diagram, and the W + W − box diagram.In all three diagrams, the virtual t quark plays a dominant role, while the c and u quarks are the secondary contributions.The FCNC process is very sensitive to the new physical effects.This suggests that it can serve as a perfect platform to search directly for new physics (NP) beyond the SM [1][2][3].
Although great progress has been made both experimentally and theoretically in the rare semileptonic decays of bottom(-strange) mesons in recent decades, those of bottomcharmed mesons have been less studied.Compared to the B (s) mesons, the B c meson is difficult to produce at the Belle experiment because the B c Bc is close to 12.5 GeV, which is far from the energy region of Υ(4S ).Moreover, according to f c / f u = (7.5 ± 1.8) × 10 −3 measured by the LHCb collaboration [109], the B c meson is also underproductivity in the pp experiment.Here, the f c and f u are the fragmentation fractions of B c and B meson, respectively, in pp collisions.As a result, the B c meson decay has received less experimental attention in the past.Recently, the LHCb collaboration reported the result of the B + c → D + s µ + µ − process [110].Using the pp collision data collected by the LHCb experiment at the center-of-mass energies of 7, 8, and 13 TeV, corresponding to a total integrated luminosity of 9 fb −1 , the LHCb collaboration did not observe significant signals in the nonresonant µ + µ − modes, but set an upper limit as f c / f u ×B(B + c → D + s µ + µ − ) < 9.6×10 −8 at the 95% confidence level.Moreover, considering that the B c → D ( * ) s ℓ + ℓ − channels have similar amounts of branching fractions [111], and the D * s needs to be reconstructed by the D s meson in the experiment, the measurement of B c → D * s ℓ + ℓ − will be more difficult.This indicates that the search for rare semileptonic decays of B c is difficult for the present experiment.However, with the high-luminosity upgrade of the Large Hadron Collider (LHC), this situation is likely to improve.In any case, with the accumulation of data in the experiment, we expect the LHCb experiment to search for these rare semileptonic decays of the B c meson.
In the theoretical sector, the rare semileptonic decays of B c have been studied by the light-front quark model (LFQM) [112], the pQCD [111], the QCD sum rule [113,114], the constituent quark model (CQM) [112], et al.The branching fractions of B c → D * s ℓ + ℓ − with ℓ = e or µ are predicted to be approximately 10 −7 .In Refs.[115][116][117], the B c → D * s µ + µ − process was been studied within the SM and beyond.In this work, we also focus on the B c → D * s ℓ + ℓ − process, where the necessary form factors are calculated via the covariant LFQM approach.To provide more physical observables, we present the angular distribution of the quasi-fourbody process B c → D * s (→ D s π)ℓ + ℓ − .The applications of the standard and(or) covariant LFQM have proved successful in the study of the meson  and baryon weak decays .The B c → D * s weak transition form factors deduced by (axial)-vector currents have been calculated in Ref. [154] with the covariant LFQM.Probably in the series of papers [124-134, 137-140, 144-153, 155, 158-181], the hadron wave function was taken as a Gaussian-like form with phenomenal parameter β, which represents the hadron structure.To fix the phenomenal parameter, the corresponding decay constant was used.However, as we all know, the decay constant is only associated with the zero-point wave function.This indicates that the oversimplified Gaussian-form wave function is not able to depict the behavior far away from the zero point.For this object, we propose to directly adopt the numerical spatial wave function by solving the Schrödinger equation with the modified Godfrey-Isgur (GI) model.By fitting the mass spectrum of the observed heavy flavor mesons, the parameters of the potential model can be fixed.This strategy avoids the β dependence, and can also reduce the corresponding uncertainty.We note that in Ref. [186], the authors used a relativistic quark model based on the quasipotential approach in QCD to study the semileptonic decay of bottom mesons.In their approach, the numerical wave functions of the mesons are obtained, thus avoiding the corresponding uncertainty.
This paper is organized as follows.After the Introduction, we illustrate the angular distributions of the quasi-four-body decays B c → D * s (→ D s π)ℓ + ℓ − (ℓ= e, µ, τ) in Sec.II.In Sec.III, we introduce the covariant LFQM and derive the formula of the weak transition form factors. Then in Sec.IV, the numerical results, including the from factors of B c → D * s and physical observables of B c → D * s (→ D s π)ℓ + ℓ − processes, are presented.Finally, this paper ends with a short summary.

II. THE ANGULAR DISTRIBUTION OF
The effective Hamiltonian associated with b → sℓ where V i j are the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements and G F = 1.16637 × 10 −5 GeV −2 [188] is the Fermi constant.Also, the C i (µ) are Wilson coefficients and the O i (µ) are four fermion operators.They all depend on the QCD renormalization scale µ.More specifically, the O c,u 1,2 are current-current operations, the O 3−6 are QCD penguin operators, the O 7,8 are electromagnetic and chromomagnetic penguin operators, and the O 9,10 are semileptonic operators, respectively.
Apart from the γ and Z penguin diagrams, and the W + W − box diagram, the long distance contribution, via the intermediate vector states (ρ, ω, ϕ.J/ψ, ψ(2S ), . . . ) (see Fig. 1) also shows an unignorable influence.By adding the factorable quark-loop contributions from O 1−6,8 to the effective Wilson coefficients C eff 7,9 , the effective Hamiltonian in Eq. (2.1) can be simplified.In the calculation, we have adopted the following effective Hamiltonian, i.e., where where the term C ′ b→sγ is the absorptive part of the b → scc → sγ rescattering [54,57,58,64,76,[189][190][191]: , and α s being adopted as α s (m b ) = 0.217 in our calculation.The short-distance contributions from the soft-gluon emission and the one-loop contributions of the four fermion operators O 1 − O 6 , and the long-distance contributions from the intermediate vector meson states are also taken into account, and have been included in the Y pert and Y res terms, respectively.The Y pert can be written as [192] Y pert ( ŝ, µ) =0.124ω( ŝ) + g( mc , ŝ)C(µ) In the Wolfenstein representation, the λ u = V ub V * us /(V tb V * ts ) can be expressed as approximately, which is a small value suppressed by λ 2 with λ = 0.22500 ± 0.00067 [188].
where m V i and Γ V i are the mass and total width of the intermediate vector meson V i respectively, and the Γ(V i → ℓ + ℓ − ) is the corresponding dilepton width.These input values are collected in Table I.In addition, the nonvanished branching fraction for the τ channel, i.e., B(ψ(2S ) → τ + τ − ) = 3.1 × 10 −3 [188], is also used.For the J/ψ and ψ(2S ) states, the small widths and the large dilepton width will have a large influence on the decay width.However, the narrow widths are also used to reject them in the experimental analysis.One the other hand, for those above the D D threshold, such as ψ(3770), ψ(4040) and ψ(4160), the board widths and mutual overlap make things difficult.Also, for the charmless vector mesons (ρ, ω and ϕ), their contributions are suppressed by the λ u factor.TABLE I: The masses, total widths and dilepton widths of the intermediate vector mesons used in Eq. (2.9).These values are quoted from the PDG [188].In this subsection, we will drive the formula of the quasifour-body decay B c → D * s (→ D s π)ℓ + ℓ − .The differential decay width of this process is where p is the four momentum of the initial B c meson, k 1 (k 2 ) and q 1 (q 2 ) are the momenta of the mesons D s (π) and the lepton ℓ − (ℓ + ), respectively, and dΦ 4 is the four-body phase space.
Taking into account the width of the D * s meson, but treating it as narrow (Γ D * s ≪ m D * s ), the width can be obtained by doing the integration as The invariant amplitude M can be calculated from where 4π , and the factor 1/2 comes from the P L(R) in the effective Hamiltonian in Eq. (2.2).
For the amplitude M D * s →D s π , it can be evaluated by the effective Lagrangian approach.The concerned effective Lagrangian is where g D * s D s π is the corresponding coupling constant.So we have the decay width of D * s → D s π as s .Obviously, the coupling constant g D * s D s π can be canceled between the vertex factor and the decay width.
Finally, with the effective Hamiltonian in Eq. ( 2.2), we can calculate the quasi-four-body decay As deduced in Ref. [196], the corresponding angular distributions can be simplified as where the explicit expressions of I i (q 2 ) and f i (θ, θ ℓ , ϕ) are shown in Table II.Compared to Ref. [196], the term I 6c is neglected since it depends on the scalar operator.As shown in Fig. 2, the θ is the angle between the −ẑ direction and pionemitted direction in the rest frame of the D * s meson, the θ ℓ is the angle made by the ℓ − with the +ẑ direction in the ℓ + ℓ − center of mass system, and the ϕ is the angle between the decay planes, i.e., the D * s → D s π plane and the virtual boson The amplitudes A L,R 0,∥,⊥ and A t are the functions of the transferred momentum square q 2 , and the seven independent form factors V, A 0,1,2 and T 1,2,3 , i.e., [57,58,196] where M ′ (M ′′ ) is the mass of the B c (D * s ) meson and mb = m b /q 2 , and (2.17) where Īi can be obtained by doing the conjugation for the weak phases of the CKM matrix elements in I i in Table II.
In addition, we should also do the following substitutions as This is the result of the operations of θ ℓ → θ ℓ − π and ϕ → −ϕ.
To separate the CP-conserving and the CP-violating effects, we define the normalized CP-averaged angular coefficients S i and the CP asymmetry angular coefficients A i as (2.20) respectively.To reduce both the experimental and theoretical uncertainties, the S i and A i have been normalized to the CP-averaged differential decay width.The other physical observables, such as the forward-backward asymmetry parameter A FB , the CP-violation A CP , and the longitudinal (transverse) polarization fractions of D * s meson F L (F T ), can thus be easily expressed in terms of these normalized angular coefficients.With the above preparations, we continue to study the physical observables.[57,58,196] in Eq. (2.15), where m2 ℓ = m 2 ℓ /q 2 and β ℓ = 1 − 4 m2 ℓ .
(b) The CP violation of the decay width can thus be estimated by (2.23) (c) The CP asymmetry lepton forward-backward asymmetry is and the CP-averaged lepton forward-backward asymmetry is The longitudinal and transverse D * s polarization fractions are Furthermore, the clean angular observables P 1,2,3 and P ′ 4,5,6,8 (more details can be found in Refs.[197,198]) are associated with the CP-averaged angular coefficients: (2.28) As pointed out in Refs.[23,57,197,198], in the large-recoiled limit, these observables are largely free of form factor uncertainties.
Finally, we also focus on the ratios, i.e., which reflect the LFU.We would like to emphasize that the lower limit of the integral of the electron mode is chosen as 4m 2 µ instead of the kinematic limit 4m 2 e in order to exclude the large enhancement dominated by the photon pole in the small q 2 region due to the C eff 7 -associated factor 1/q 2 .In the B → Kℓ + ℓ − process, the experimental measurements of the ratio R eµ K by Belle [9][10][11][12] and BABAR [15] are in agreement with the SM prediction, while the LHCb result [19,21,22] shows a clear deviation from the SM expectation (see Fig. 4 of Ref. [22]) with 3.1σ.We note that in Ref. [199], the authors used the ratios R τµ K ( * ) to study the LFU violation, and found that they can deviate from the SM prediction even if the NP couplings are universal.Therefore, in order to use these ratios to study the LFU violation, we should compare the allowed ranges, considering both the solutions with only universal couplings and those with universal and nonuniversal components.Whatever, the ratio in the B c → D * s ℓ + ℓ − sector is also interesting to investigate whether it is consistent with the SM expectation or not.The breaking of the LFU may require an expansion of the gauge structure of the SM, and of course probes the NP effects [200].

III. WEAK TRANSITION FORM FACTORS
The standard and(or) covariant LFQMs have been widely used to study the decays of mesons  and baryons .In the conventional LFQM framework, the consistent quark (or antiquark) of the meson is required to be on its mass shell, and thus the initial (or final) meson is offshell.This procedure misses the zero-mode effects and makes the matrix element noncovariant.To avoid this shortcoming, Jaus [118,122] proposed a covariant framework for the S -waved pseudoscalar and vector meson decays in which the zero-mode contributions are systematically taken into account.Cheng et al. [124,132] extended this approach to the case of the P-wave meson (such as scalar, axial-vector and tensor mesons).The physical quantities, such as the decay constant and the form factor of the weak transition, are obtained in terms of the Feynman loop integration.Unlike the conventional LFQM, the covariant LFQM requires the initial (or final) meson to be on its mass shell.For more details on the difference, see Refs.[124,146].In this section, we will use the covariant LFQM to calculate the B c → D * s form factors.Following Ref. [154], the B c → D * s weak transition form factors deduced by (axial-)vector currents are defined as where we use the convention ϵ 0123 = +1 and define P µ = p ′ µ + p ′′ µ and q µ = p ′ µ − p ′′ µ , and ε is the polarization vector of the D * s meson.These amplitudes can also be parametrized as the Bauer-Stech-Wirbel (BSW) form [201], i.e., with M ′ (M ′′ ) being the mass of the parent (daughter) meson.These two definitions are related by the relations [154] In addition, the (pseudo)tensor current amplitudes can be defined as [202,203] where, we have T 1 (0) = T 2 (0) since the identity 2σ µν γ 5 = −iϵ µναβ σ αβ .
The form factors require a nonperturbative calculation.In this work, we use the covariant LFQM to calculate the relevant form factors for the weak transition.In this approach, the constituent quark and the antiquark inside a meson are off shell.We define the incoming (outgoing) meson to have the momentum , where p ′(′′) 1 and p 2 are the off-shell momenta of the quark and the antiquark, respectively.These momenta can be expressed in terms of the internal variables (x i , ⃗ k ′ ⊥ ) (i = 1, 2), defined by (3.5)They must also satisfy x 1 + x 2 = 1.
The one-loop Feynman diagram for meson weak transition amplitude, where P ′ (P ′′ ) is the momentum of the incoming (outgoing) meson, p ′(′′) 1 and p 2 are the momenta of the quark and antiquark, respectively.The symbol "cross" denotes the weak interaction vertex.
According to Refs.[124,154], the corresponding weak transition matrix element at the one-loop level can be calculated in terms of the Feynman loop integral, as shown in Fig. 3. Then the form factors can be extracted from the corresponding matrix element.To write down the transition amplitude, we need the meson-quark antiquark vertices for the initial meson as iΓ ′ = H ′ P γ 5 , and that of the outgoing meson as i(γ 0 Γ ′′ † γ 0 ) with [124,154], where the subscripts P and V denote the pseudoscalar and vector meson, respectively.
For Fig. 3, the concrete expression of the transition amplitude for P → V can be expressed as where N ′(′′) and N 2 = p 2 2 − m 2 2 come from the propagators of the quarks.The superscripts V, A, T , and T 5 represent the vector, axial-vector, tensor, and pseudotensor currents, respectively.The traces S V µν are written as (3.7) To make reading easier, the relevant expressions of the traces S A,T,T 5 µν are collected in Appendix A.
Following Refs.[123,124,132], the execution of the p ′− 1 integration went to the replacement: where we define with To write down the concrete expression of Ŝ µν , we should take into account the so-called zero-mode contribution.As shown in Refs.[124,141], after doing the integration in Eq. (3.8) we have p 2 = p2 , and where ω = (2, 0, ⃗ 0 ⊥ ) is a lightlike four vector in the lightfront coordinate.Following the discussions in a series of papers [123,124,132,141], for avoiding the ω dependence, we need to do the following replacements [123,124,132]: in Eqs.(3.7), (A.1), (A.4), and (A.5).Here, (3.12) After performing the replacements (3.11) in the decay amplitudes (3.7) and (A.1), the form factors g, f , a + , and a − can be obtained from the terms proportional to the ϵ µναβ P α q β , g µν , P µ P ν and P µ q ν , and q µ P ν and q µ q ν , respectively.The ϵ * µ V P ′′ µ = 0 is used here.Finally, the expressions of these form factors in covariant LFQMs can be written as [122,124,154] The form factors deduced by (axial)vector currents defined in Eq. (3.2) can thus be evaluated by (3.17) Analogously, we can obtain the concrete expressions of the (pseudo)tensor form factors defined in Eq. (3.4) as [132] T ] Following the treatment in Ref. [124], h M is taken as where M′(′′) , and ϕ s is the space wave function of the pseudoscalar or vector meson.
In the previous theoretical work [124,154], the phenomenological Gaussian-type wave functions are widely used.It inevitably introduces the dependence of the parameter β.The phenomenological parameter β can be fixed by the decay constant [123,124,132].However, as we all know, the decay constant is only associated with the meson wave function at the end point q 2 = 0.This indicates that the simple wave function Eq. (3.22) deviating from the q 2 = 0 region may be unreliable.Taking advantage of the modified GI model [204], we can obtain the numerical spatial wave functions of the mesons concerned.By replacing the form in Eq. (3.22) with where c n are the expansion coefficients of the corresponding eigenvectors and l is the orbital angular momentum of the meson, we can avoid the corresponding uncertainty.In Table III, we collect the expansion coefficients c n of the meson wave functions involved.In addition, the factor √ 4π is needed to satisfy the normalization: Besides, the R nl is the simple harmonic oscillator wave function as (3.26)The parameter β = 0.5 GeV in the above equation is consistent with Ref. [204].

A. The form factors
With the input of the numerical wave functions, and the concrete expressions of the seven form factors in Eqs.(3.13)- (3.20), we present in this subsection the numerical results of B c → D * s form factors.Following the approach described in Refs.[122,124], we assume the condition q + = 0.This implies that our form factor calculations are performed in the spacelike region (q 2 < 0), and therefore we need to extrapolate them to the timelike region (q 2 > 0).To perform the analytical continuation, we utilize the z-series parametrization [141] where a i (i = 1, 2) are free parameters needed to fit in the q 2 < 0 region, and the z(q 2 ) is taken as The calculated masses and the expansion coefficients c n of the wave function of the mesons involved [204].The masses are given in units of MeV.

B. The angular distributions and physical observables
With the above preparations, in this subsection we present our numerical results of the branching fractions and some angular observables, i.e., the CP-averaged normalized angular coefficients S i , the lepton's forward-backward asymmetry parameter A FB , and the longitudinal (transverse) polarization fractions of the D * s meson F L(T ) .In addition, we also investigate the clean angular observables P 1,2,3 and P ′ 4,5,6,8 .The hadron and lepton masses are quoted from the PDG [188], as well as the lifetime τ B c = 0.510ps and the branching fraction B(D * s → D s π) = 5%.First, we focus on the angular coefficients S i and A i defined in Eq. (2.20).The q 2 dependence of the normalized CP-averaged angular coefficients S i are presented in Fig. 5, while the CP asymmetry angular coefficients A i are shown in Fig. 6.The blue dashed lines and the magenta solid lines represent the muon and the tau channels, respectively.Since the electron channel shows similar behavior to the muon channel, we will only present our results for the muon and the tau channels here.In the energy regions of 8.0 < q 2 < 11.0 and 12.5 < q 2 < 15.0 GeV 2 , we use the gray areas to mark the contributions from the charmonium states J/ψ and ψ(2S ).In our calculation, we adopted phenomenological and model-dependent treatment, i.e., the Breit-Wigner ansatz to model the corresponding contribution.In the experiment, these two regions are generally truncated.The CP asymmetry angular coefficients, A i , are shown to be very small in the SM, due to the direct CP violation being proportional to the Im[V ub V * us /V tb V * ts ], which is around 10 −2 .This character is very clear in Fig. 6.Also, the S 7,8,9 are also very small compared to the other angular coefficients S i .These angular coefficients are important physical observables to reveal the underlying decay mechanism, and can be checked by future measurements at LHCb.
We further evaluate the CP-averaged differential branching fractions by using Eqs.(2.21) and (2.22).The q 2 dependence of the differential branching fractions are shown in Fig. 7, where the red, blue and magenta curves represent the e, µ, and τ modes, respectively.The gray areas also denote the charm loop contributions from the charmonium states J/ψ and ψ(2S ).In Table VI, we present our result of the branching fractions and their ratios in different q 2 bins.In the four q 2 intervals, i.e., [1.1, 6.
muon modes can reach up to 10 −8 , and the ratio R eµ = 1, which is consistent with the SM prediction and reflects the LFU.In the high q 2 region, that is [15.0, 17.0] GeV 2 , the branching fraction of the tau mode is on the order of magnitude of 10 −9 .We also obtain the ratio R τµ = 0.384.In the region of 1.1 < q 2 < 6.0 GeV 2 , we have the branching fractions as In addition, combined with the branching fraction B(D * s → D s π) = 5%, we have which may well be tested by the ongoing LHCb experiment.We also investigate the physical observables, i.e., the lepton forward-backward asymmetry parameter A FB and the longitudinal (transverse) polarization fractions F L (F T ).The q 2 TABLE VI: Our results of the branching fractions of B c → D * s (→ D s π)ℓ + ℓ − (ℓ= e, µ, τ) (in units of 10 −8 ) in different q 2 bins.dependence of these observables is presented in Figs. 8 and  9, respectively.Their averaged values in different q 2 bins, defined by FIG. 5: The q 2 dependence of normalized CP-averaged angular coefficients S i , where the blue dashed and magenta solid curves are our results for the µ and τ modes, respectively.
FIG. 6: The q 2 dependence of normalized CP asymmetry angular coefficients A i , where the blue dashed and magenta solid curves are our results for the µ and τ modes, respectively.where A = (A FB , F L , F T ), are shown in Table VII.In addition, we present our results for the q 2 dependent clean angular observables P 1,2,3 and P ′ 4,5,6,8 in Fig. 10.In Ref. [23], the LHCb collaboration reported the measurement of the form-factor-independent observables P ′ 4,5,6,8 of the B 0 → K * 0 µ + µ − decay.In particular, in the interval of 4.30 < q 2 < 8.68 GeV 2 , the observable P ′ 5 shows 3.7σ discrepancy with the SM prediction [24].After integration over the energy region 1.0 < q 2 < 6.0 GeV 2 , the discrepancy is determined to be 2.5σ.So we want to investigate these clean angular observables in the rare semileptonic decay of bottomcharmed meson.In order to exclude the charmonium contributions and make it easy to check experimentally, we also present the averaged values of these observables in different q 2 intervals in Table VIII.The averaged value in a q 2 bin is defined by Eq. (4.3).
In general, this quasi-four-body decay provides a set of physical observables to study the corresponding weak interaction, and in particular the ratios of the branching fractions R eµ and R τµ , as well as the clean angular coefficients P i and P ′ j , can be helpful to search for the NP effects beyond the SM.We call for the ongoing LHCb experiment to search for this process and to measure the corresponding physical observables.

V. SUMMARY
In this work, we have studied the B c → D * s transition form factors deduced by the (axial)vector and (pseudo)tensor currents, and, in the future, investigate the angular distributions of the quasi-four-body processes To describe the weak process, the relevant seven independent form factors are calculated by utilizing the covariant LFQM approach.The concerned meson wave functions are adopted as the numerical wave functions, which are extracted from the solution of the modified GI model.This treatment avoids the β dependence and thus reduces the corresponding uncertainty.Our results of form factors are compared with other approaches.In particular, for the (pseudo)tensor currents deduced form factors T 1,2,3 (q 2 = 0), our results agree with the pQCD prediction.More theoretical works, especially the LQCD and QCD sum rule (or light-cone sum rule) calculation, are highly appreciated to test our result and to refine the corresponding topic.
With the obtained form factors, the rare semileptonic decays B c → D * s (→ D s π)ℓ + ℓ − are studied.Not only the branching fractions, the lepton-side forward-backward asymmetry parameter A FB , and the longitudinal and transverse polarization fractions F L and F T , but also the angular coefficients S i and A i are investigated.Numerically, the concerned cascade decays with e or µ final states are around 10 −8 , which need to be tested by other approaches and ongoing experiments.Moreover, the A FB and F L(T ) are important physical observables, and they are also feasible observables in the future LHCb experiment, so we look forward to the experimental results.In addition, the ratios of the branching fractions are also calculated to validate whether or not the LFU violated.Furthermore, the clean coefficient observables P i and P ′ j are presented, which reduce the uncertainty from the form factors and can be a possible signal to search for the NP effects.Since these observations are largely free of form factor uncertainties in the large-recoiled limit, and are feasible to measure experimentally, we strongly encourage our experimental colleagues  to measure them.
Overall, in this work we have systematically studied the angular distribution of B c → D * s (→ D s π)ℓ + ℓ − (ℓ= e, µ, τ) with the form factors obtained by the covariant LFQM.We live in the hope that with the completion of the LHCb experiment prepared for the run 3 and run 4 of the LHC and the improvement of the experimental capabilities, this rare semileptonic process can be discovered and we expect that the predicted physical observables can be tested.
0], [6.0, 8.0],[11.0,12.5], and [15.0, 17.0] GeV 2 , the branching fractions of the electron and FIG. 4: The q 2 dependence of the B c → D * s weak transition form factors. Here, the four dependent form factors deduced by (axial)vector are presented in the left panel, while the three dependent ones deduced by (pseudo)tensor are shown in the right panel.

FIG. 8 :
FIG.8:The q 2 dependence of lepton forward-backward asymmetry parameter A FB in B c → D * s (→ D s π)ℓ + ℓ − [ℓ=e (left panel), µ (center panel), and τ (right panel)] processes, where the red, the blue, and the magenta curves are our results from the e, µ, and τ modes, respectively.

FIG. 9 :
FIG.9:The q 2 dependence of D * s longitudinal (transverse) polarization fractionsF L (F T ) in B c → D * s (→ D s π)ℓ ℓ − [ℓ=e (left panel), µ (center panel), and τ (right panel)] processes, where the red, the blue, and the magenta curves are our results from the e, µ, and τ modes, respectively, and the solid and dashed curves represent the F L and F T , respectively.

TABLE II :
The explicit expressions of the angular coefficients I i and f i

TABLE IV :
Our results of the weak transition form factors of B c → D * s by using the covariant LFQM.

TABLE V :
Theoretical predictions of the B c → D * s transition form factors at the end point q 2 = 0 using different approaches.

TABLE VII :
The averaged forward-backward asymmetry ⟨A FB ⟩ and the longitudinal (transverse) polarization fractions ⟨F L ⟩(⟨F T ⟩) in different q 2 bins.