Probing the $P$-wave charmonium decays of $B_c$ meson
Zhou Rui
hep-ph hep-ex









jindui1127@126.com
^1College of Sciences, North China University of Science and Technology,
                          Tangshan 063009,  China





Motivated by the large number of B_c meson decay modes observed recently by  several detectors at the LHC,
we present a detailed analysis of the B_c meson decaying to the P-wave charmonium states and a light pseudoscalar
 (P) or vector (V) meson within the framework of perturbative QCD factorization.
 The P-wave charmonium distribution amplitudes are extracted from the n=2,l=1 Schrödinger states for a Coulomb potential,
 which can be taken  as the universal nonperturbative objects to analyze the hard exclusive processes with P-wave charmonium production.
It is found that these decays have  large branching ratios of the order of 10^-5∼ 10^-2,
 which seem to be in the reach of future experiments.
  We also provide predictions for the polarization fractions and relative phases  of B_c→ (χ_c1,χ_c2,h_c)V decays.
It is expected  that  the longitudinal polarization amplitudes  dominate the branching ratios according to the quark helicity analysis,
and the magnitudes and phases of parallel polarization amplitude are approximately equal to the perpendicular ones.
The obtained results are compared with available experimental data, our previous studies, and numbers from other approaches.


13.25.Hw, 12.38.Bx, 14.40.Nd 




Probing  the P-wave charmonium decays of B_c meson
    Zhou Rui^1
    December 30, 2023
==================================================






§ INTRODUCTION




In the quark model, P-wave charmonium states are expected as the orbital excitation of the cc̅  assignments
with the orbital angular momentum L=1.
 Since the charm-anticharm quarks pair can be in the spin singlet or spin triplet states,
in terms of the spectroscopic notation ^2s+1L_J, there are four types of P-wave charmonium states, namely, χ_c0(^3P_0),
χ_c1(^3P_1), χ_c2(^3P_2), and h_c(^1P_1).
The current experimental knowledge of these P-wave charmonium states is summarized in Table <ref> <cit.>.
 Experimentally the productions of these P-wave charmonium states have been seen in the   hadronic B decays:
 B→χ_c1π <cit.>,
B→χ_c0 K^* <cit.>, B→χ_c1,c2 K^(*) <cit.>,
B→ h_c K^(*) <cit.>,  B_s→χ_c1ϕ <cit.>, 
B→χ_c1,c2Kπ <cit.>, B→χ_c1,c2ππ K <cit.>,
and  in Λ_b^0 baryon decay: Λ_b^0→χ_c1,c2pK^- <cit.>.
As for   hadronic B_c decays,  the first evidence of B^+_c→χ_c0(→ K^+K^-)π^+  <cit.>
is reported with a significance of 4.0 standard deviations by  the LHCb experiment.
The measured product of the ratio of cross sections and branching fraction is

    σ_B_c^+/σ_B^+×ℬ(B^+_c→χ_c0π^+)
    =(9.8^+3.4_-3.0(stat)± 0.8(syst))× 10^-6,

where  σ_B_c^+(σ_B^+) is the production cross sections for B_c(B) meson.

As is well known, the pseudoscalar B_c is composed of two  heavy-flavored quarks
and thus represents a unique laboratory  to study heavy-quark dynamics and CP violation.
Since each of the two heavy quarks can decay with the other as a spectator,
the B_c meson has rich decay channels, and offers a promising  opportunity
to study nonleptonic weak decays of heavy mesons, to test the standard model (SM), and even to reveal any new physics beyond SM.
Decays of B_c mesons to the final states including a
charmonium meson are of special interest. First, these decay
modes   provide a sensitive laboratory for studying strong
interaction   effects in a heavy meson system.
Second, those decays involve two energy scales, the bottom quark mass m_b and charm quark mass m_c.
The higher order corrections within the framework of quantum chromodynamics (QCD), described by the expansion  of m_c/m_b rather than Λ_QCD/m_b with Λ_QCD is the QCD scale, may be relatively large, and therefore are more subtle in theoretical studies.
 Third, one can search for charmonium and charmoniumlike exotic states in one of
the intermediate final states such as χ_cJπ and χ_cJππ (J=0,1,2), which may be  important to understand the detailed
dynamics of the multibody B_c decay modes.
Besides, such B_c decays provide a direct probe of charmonium properties
by reconstructing the charmonium state from its decay to some known final state.


Phenomenologically the B_c meson  decays  into  various charmonium states  have been widely studied
in the literature. Earlier, a lot of work has been done in the semileptonic
and nonleptonic <cit.> decays of the meson B_c to S-wave charmonium
mesons. Also, the P-wave charmonium  decays of B_c meson
have been considered previously by other authors
<cit.>.
Furthermore, some semileptonic and nonleptonic decays of B_c  into the D-wave charmonium states
have been analyzed  in the framework of the instantaneous Bethe-Salpeter method <cit.>.
More recently, the exclusive decays of the B_c meson into P-wave orbitally excited charmonium and
a light meson have been investigated using the nonrelativistic QCD effective theory <cit.>,
 where the next-to-leading order relativistic corrections to the corresponding form factors are considered.


As a successive work of <cit.>,
in the present work we will focus on the B_c decays
involving a P-wave charmonium state  and a light pseudoscalar or vector
meson  in the final states employing the Perturbative QCD (PQCD) approach based on the k_T factorization theorem.
Similar to the case of S-wave charmonium states <cit.>,
the P-wave charmonium  distribution amplitudes (DAs)
  can also be expressed  as  an   associated factor,  extracted from the P-wave Schrödinger states for a Coulomb potential,
   multiply by the asymptotic models of the corresponding twists for light mesons.
   With the help of  P-wave DAs, we can make quantitative predictions here, and
 provide a ready reference to existing and forthcoming experiments.

The rest of this article is organised in the following way.
In Sec. <ref>, the Hamiltonian and kinematics, and
the P-wave charmonium DAs  are shown in cases of   scalar, axial-vector, and tensor states. Then
the calculations of these decay amplitudes in the PQCD framework are briefly reviewed.
In Sec. <ref>, the adopted parameters, numerical results and discussion are given in detail.
Finally, the conclusions are given in Sec. <ref>.
The evaluation of the P-wave charmonium distribution
amplitudes is relegated to the Appendix.



§ FORMALISM 

 

 §.§ HAMILTONIAN AND KINEMATICS 


 The effective Hamiltonian describing the B_c nonleptonic decays into charmonium and a light pseudoscalar or vector meson
is given by <cit.>
 
    ℋ_eff=G_F/√(2)V^*_cbV_uq[C_1(μ)O_1(μ)+C_2(μ)O_2(μ)],

 where q=s,d stands for a down type light quark. G_F is the Fermi constant. V^*_cb and V_uq are the Cabibbo-Kobayashi-Maskawa
(CKM) matrix elements.
 C_1,2(μ) are the perturbatively calculable Wilson coefficients,
 which encode the short-distance effects above the renormalization scale μ,
 while  O_1,2(μ) are the corresponding  local four-quark operators,
whose expressions  read as
 
    O_1(μ)   =   b̅_αγ^ν(1-γ_5)c_β⊗u̅_βγ_ν(1-γ_5)q_α,
    
    O_2(μ)   =   b̅_αγ^ν(1-γ_5)c_α⊗u̅_βγ_ν(1-γ_5)q_β,

where α and β are color indices and the summation convention
over repeated indices is understood.
Since the  Hamiltonian involves  four different flavor quarks, it means that these decays are uncontaminated by
the contributions from the penguin operators, and thus the direct CP asymmetries are absent naturally.

The calculation is carried out in the rest frame of B_c meson, the  B_c meson momentum P_1,
the recoiled charmonium meson momentum P_2, and the ejected light meson momentum P_3 
are defined in the light cone coordinates as

    P_1   =   M/√(2)(1,1,0_ T),   P_2=M/√(2)(1-r_3^2,r_2^2,0_ T),    P_3=M/√(2)(r_3^2,1-r_2^2,0_ T),

with the mass ratio r_2,3=m_2,3/M and M (m_2) is  the B_c (charmonium) meson mass, while
m_3 is the (chiral) mass of the  (pseudoscalar) vector meson.
The momentum of the valence quarks k_1,2,3, whose notation is displayed in Fig. <ref>, is parametrized as

    k_1   =   x_1P_1+k_ 1T,   k_2=x_2P_2+k_ 2T,    k_3=x_3P_3+k_ 3T,

where k_iT , x_i represent the transverse momentum and longitudinal momentum fraction of the quark/anti-quark inside the meson.
When the final states contain a axial-vector charmonium  and a vector meson, the longitudinal
polarization vectors ϵ_L and transverse polarization vectors ϵ_T can be defined as

    ϵ_2L   =   1/√(2(1-r_3^2))r_2(1-r_3^2,-r_2^2,0_ T),   ϵ_2T=(0,0,1_ T),
    ϵ_3L   =   1/√(2(1-r_2^2))r_3(-r_3^2,1-r_2^2,0_ T),   ϵ_3T=(0,0,1_ T),

which satisfy the normalization ϵ^2_L= ϵ^2_T=-1
and the orthogonality ϵ_2L· P_2=ϵ_3L· P_3=0.

For a tensor charmonium, the polarization tensor ϵ_μν(λ) with
helicity λ can be constructed via the polarization vector ϵ_μ  <cit.>:

    ϵ_μν(± 2)   =   ϵ_μ(±)ϵ_ν(±), 
    ϵ_μν(± 1)   =   1/√(2)[ϵ_μ(±)ϵ_ν(0)+ϵ_ν(±)ϵ_μ(0)], 
    ϵ_μν(0)   =   1/√(6)[ϵ_μ(+)ϵ_ν(-)+ϵ_μ(-)ϵ_ν(+)]+
    √(2/3)ϵ_μ(0)ϵ_ν(0),

with ϵ(±)=ϵ_2T and ϵ(0)=ϵ_2L.
It is convenient to define another polarization vector ϵ_∙μ(λ)=
m_2ϵ_μν(λ)v^ν/P_2· v, which satisfy

    ϵ_∙μ(± 2)=0,  ϵ_∙μ(± 1)= m_2√(1/2)ϵ_2L· v/P_2 · vϵ_2Tμ,  ϵ_∙μ(0)= m_2√(2/3)ϵ_2L· v/P_2 · vϵ_2Lμ.

The contraction is evaluated as ϵ_2L· v/P_2· v=1/m_2 by
neglecting the light meson mass, then we get the relations
ϵ_∙ T  =ϵ_∙(±1)= √(1/2)ϵ_2T and
ϵ_∙ L  =ϵ_∙(0)=√(2/3)ϵ_2L.
Note that ϵ_∙ has the same energy scaling as the
usual polarization vector of a vector meson. It makes the calculations of B_c decays into a tensor meson are similar to
those of vector analogues by replacing the  polarization vector with the corresponding ϵ_∙.
 

 




 §.§  MESONS WAVE FUNCTION AND THE DISTRIBUTION AMPLITUDES 

In the considered decays, there are three typical scales: M, m_2, and the heavy-meson and heavy-quark mass
difference Λ̅. They allow for a consistent power expansion in m_2/M and in Λ̅/m_2 under the hierarchy of
Λ̅≪ m_2≪ M. In the heavy-quark and large-recoil limits, based on the k_T factorization theorem,
the decay amplitudes are expressed as the convolution of the hard kernels  with the relevant meson wave functions.
The hard kernels can be treated by perturbative QCD at the leading order in an α_s expansion (single gluon exchange
as depicted in Fig. <ref>). The higher-order radiative corrections generate the logarithm divergences, which can be absorbed
into the meson wave functions. One also encounters double logarithm divergences when collinear and soft divergences
overlap, which can be summed to all orders to give a Sudakov factor. After absorbing all the soft dynamics, the initial and
final state meson wave functions can be treated as nonperturbative inputs, which are not calculable but universal.

Analogous  to the B meson <cit.>,
up to first order in 1/M under above hierarchy,
the B_c meson wave functions are decomposed into the following Lorentz structures <cit.>:

    ∫ d^4z e^ik_1· z⟨ 0|b̅_α(0)c_β(z)|B_c(P_1)⟩=
    i/√(2N_c){(P_1+M)γ_5[ϕ_B_c(k_1)-/n-/v/√(2)ϕ̅_B_c(k_1)]}_βα

with the two lightlike vectors n=(1,0,0_ T) and v=(0,1,0_ T).
N_c=3 is the color factor.
Here, we only consider the  dominant Lorentz structure from the first term, while the second Lorentz structure
 starting from the next-to-leading-power  Λ̅/M is numerically neglected <cit.>.
 In coordinate space the distribution amplitude ϕ_B_c is adopted in the form <cit.>
 
    ϕ_B_c(x)=N_Bx(1-x)e^-m_b+m_c/8m_bm_cω(m_c^2/x+m_b^2/1-x),

with the shape parameter ω=0.5 ± 0.1 GeV   related to  the  factor N_B   by the normalization
 
    ∫^1_0ϕ_B_c(x)dx=1.


For the P-wave charmonium states, we use the abbreviations A, S, and T correspond to axial-vector, scalar,
and tensor charmonium meson, respectively. In terms of
the notation in Ref <cit.>, the nonlocal matrix element for the  longitudinally and transversely polarized
axial-vector  and scalar charmonium meson can be decomposed as

    ⟨ A (P_2,ϵ_2L)|c̅_α(z)c_β(0)|0⟩   =   1/√(2N_c)∫_0^1dx e^ixP_2· z
    [m_2γ_5/ϵ_2Lψ^L_A(x)+γ_5/ϵ_2L/P_2 ψ_A^t(x)]_βα,
    ⟨ A (P_2,ϵ_2T)|c̅_α(z)c_β(0)|0⟩   =   1/√(2N_c)∫_0^1dx e^ixP_2· z
    [m_2γ_5/ϵ_2Tψ_A^V(x)+γ_5/ϵ_2T/P_2ψ_A^T(x)]_βα,
    ⟨ S (P_2)|c̅_α(z)c_β(0)|0⟩   =   1/√(2N_c)∫_0^1dx e^ixP_2· z
    [/P_2ψ_S^v(x)+m_2ψ_S^s(x)]_βα,

where the DAs ψ^L,T,v(x) are of twist-2, and ψ^t,V,s(x) of twist-3.
As mentioned  in the Introduction,
the charmonium DAs are parameterized using a combination of an universal  factor 𝒯(x) for the P-wave states
 and the asymptotic models Φ_asy(x), given by

    ψ^i∝Φ_asy^i(x)𝒯(x).

The expression of 𝒯(x) can be extracted from
P-wave Schrödinger states for a Coulomb potential, which are derived
in the Appendix. The asymptotic forms of Φ_asy^i(x) in Eq. (<ref>)
for the axial-vector mesons can be related to the ones calculated in QCD sum rules by <cit.>

    Φ_A^L(x)=f_A/2√(2N_c)ϕ_∥(x),   Φ_A^T(x)= f^⊥_A/2√(2N_c)ϕ_⊥(x),  Φ_A^t(x)=f^⊥_A/2√(2N_c) h^(t)_∥(x),   Φ_A^V(x)=f_A/2√(2N_c) g^(a)_⊥(x),

where f_A (f_A^⊥ ) is the vector  (tensor) decay constants.
The leading twist ϕ_∥,⊥ can be expanded in a series of Gegenbauer polynomials <cit.>

    ϕ_∥   =   6x(1-x)[a_0^∥+3a_1^∥(2x-1)+⋯],
    ϕ_⊥   =   6x(1-x)[a_0^⊥+3a_1^⊥(2x-1)+⋯].

Owing to the G-parity, ϕ_∥ (ϕ_⊥) for ^3P_1  state  is symmetric
(antisymmetric) under the exchange of quark and antiquark momentum fractions in the SU(3) limit.
On the contrary, ϕ_∥ is antisymmetric for ^1P_1 states, while ϕ_⊥
is symmetric in this case.
Thus the  asymptotic forms for twist-2 can be written as
 
    ϕ_∥(x)   =   6a_0^∥x(1-x),   ϕ_⊥(x)=18a_1^⊥x(1-x)(2x-1)     for  ^3P_1; 
    ϕ_∥(x)   =   18a_1^∥x(1-x)(2x-1),   ϕ_⊥(x)=6a_0^⊥x(1-x)    for  ^1P_1.

Neglecting the three-parton distribution amplitudes containing gluons and terms proportional to
light quark masses,
the twist-3 DAs can be related to the twist-2 ones by Wandzura-Wilczek-type relations <cit.>:

    g_⊥^(a)(x)   =   1/2[∫_0^x du ϕ_∥(u)/1-u+∫_x^1 dv ϕ_∥(u)/u], 
    
    h_∥^(t)(x)   =   (2x-1)[∫_0^x du ϕ_⊥(u)/1-u-∫_x^1 dv ϕ_⊥(u)/u],


which further give
 
    h_∥^(t)(x)   =   3a_1^⊥(2x-1)(1-6x+6x^2),    g_⊥^(a)(x)=3/4a_0^∥(1+(2x-1)^2)     for  ^3P_1; 
    
    h_∥^(t)(x)   =   3a_0^⊥(2x-1)^2,             g_⊥^(a)(x)=3/2a_1^∥(2x-1)^3        for  ^1P_1.

Combining  Eqs. (<ref>), (<ref>), and  (<ref>), we derive

    ψ^L(x)   =   f_A/2√(2N_c)N_Lx(1-x)𝒯(x),
    ψ^T(x)   =   f_A^⊥/2√(2N_c)N_Tx(1-x)(2x-1)𝒯(x),
    ψ^t(x)   =   f_A^⊥/2√(2N_c)N_T/6(2x-1)[1-6x+6x^2]𝒯(x),
    ψ^V(x)   =   f_A/2√(2N_c)N_L/8[1+(1-2x)^2]𝒯(x),

for ^3P_1 states, and

    ψ^L(x)   =   f_A/2√(2N_c)N_Tx(1-x)(2x-1)𝒯(x),
    ψ^T(x)   =   f_A^⊥/2√(2N_c)N_Lx(1-x)𝒯(x),
    ψ^t(x)   =   f_A^⊥/2√(2N_c)N_L/2(1-2x)^2𝒯(x),
    ψ^V(x)   =   f_A/2√(2N_c)N_T/12(2x-1)^3𝒯(x),

for  ^1P_1 states.
Note that the Gegenbauer moments a_0,1^∥,⊥ are absorbed into the coefficients N_L,T
to satisfy the  normalization conditions <cit.>

    ∫_0^1N_Lx(1-x)𝒯(x)dx=1,   ∫_0^1N_Tx(1-x)(2x-1)^2𝒯(x)dx=1.


Because of the charge conjugation invariance, twist-2 and twist-3 DAs of the scalar meson should satisfy
 ψ_S^v(x)=-ψ_S^v(1-x) and ψ_S^s(x)=ψ_S^s(1-x), respectively <cit.>.
In general, the  asymptotic twist-2 DA can also be expanded in the Gegenbauer polynomials with only odd component contribute <cit.>. Based on the description of the charmonium states DAs given above,
ψ^v_S(x) can be recast to the form

    ψ^v_S(x)   =   f_S/2√(2N_c)N_Tx(1-x)(2x-1)𝒯(x).

As for the  twist-3 DA, we adopt the same  asymptotic form as the pseudoscalar mesons <cit.>:

    ψ^s_S(x)   =   f_S/2√(2N_c)N_S𝒯(x),

with the  normalization condition

    ∫_0^1ψ^s_S(x)dx   =   f_S/2√(2N_c).


The nonlocal matrix element associating with the tensor charmonium can be decomposed as <cit.>

    ⟨ T (P_2,ϵ_∙ L)|c̅_α(z)c_β(0)|0⟩   =   1/√(2N_c)∫_0^1dx e^ixP_2· z
    [m_2/ϵ_∙ Lψ_T(x)+/ϵ_∙ L/P_2 ψ^t_T(x)]_βα,
    ⟨ T (P_2,ϵ_∙ T)|c̅_α(z)c_β(0)|0⟩   =   1/√(2N_c)∫_0^1dx e^ixP_2· z
    [m_2/ϵ_∙ Tψ^V_T(x)+/ϵ_∙ T/P_2 ψ^T_T(x)]_βα,

for the longitudinal    and transverse polarizations, respectively.
ψ_T(x), ψ^T_T(x)  are leading twist   DAs, and ψ^V_T(x), ψ^t_T(x) are twist-3 ones,
which are related to the ones given in Ref. <cit.>

    Φ_T(x)=f_T/2√(2N_c)ϕ_∥(x),   Φ_T^T(x)= f^⊥_T/2√(2N_c)ϕ_⊥(x),  Φ_T^t(x)=f^⊥_T/2√(2N_c) h^(t)_∥(x),   Φ_T^V(x)=f_T/2√(2N_c) g^(v)_⊥(x).

In SU(3) limit, due to the G-parity of the tensor meson,
all of the DAs are antisymmetric under the replacement x→ 1-x.
Here we take the following approximate forms of twist-2 as <cit.>


    ϕ_∥(x)=ϕ_⊥(x)=N_Tx(1-x)(2x-1),

and the corresponding expressions for the twist-3 DAs can be derived through the Wandzura-Wilczek relations  as <cit.>

    h^(t)_∥(x)=N_T/4(2x-1)(1-6x+6x^2),   g^(v)_⊥(x)=N_T/6(2x-1)^3.

Now we can collect the DAs of tensor charmonium states below:

    ψ_T(x)   =   f_T/2√(2N_c)N_Tx(1-x)(2x-1)𝒯(x),
    ψ^T_T(x)   =   f_T^⊥/2√(2N_c)N_Tx(1-x)(2x-1)𝒯(x),
    ψ^t_T(x)   =   f_T^⊥/2√(2N_c)N_T/4(2x-1)[1-6x+6x^2]𝒯(x),
    ψ^V_T(x)   =   f_T/2√(2N_c)N_T/6(2x-1)^3𝒯(x),

with the normalization conditions <cit.>

    ∫_0^1(2x-1)ψ^(T)_T(x)=f_T^(⊥)/2√(2N_c).


For the wave functions of light pseudoscalar and vector mesons,
the same forms and parameters are adopted as <cit.> and one
is referred to the original literature <cit.>.



 §.§  THE DECAY AMPLITUDES 

In the PQCD approach, the decay amplitudes are expressed
as the convolution of the hard kernels H with the relevant
meson wave functions Φ_i

    𝒜(B_c→ M_2M_3) = ∫ d^4k_1d^4k_2d^4k_3 Tr[C(t)Φ_1(k_1)Φ_2(k_2)Φ_3(k_3)
    H(k_1,k_2,k_3,t)].

“Tr” denotes the trace over all Dirac structures and color indices.
The hadron wave functions Φ_i absorbed all the nonperturbative components have been described in  Sec. <ref>.
The hard kernel H(k_1,k_2,k_3,t) describes the four quark
operator and the spectator quark connected by a hard gluon,
which can be perturbatively calculated including all possible
Feynman diagrams without end-point singularity.
In the following, we start to compute the decay amplitudes of the concerned  decays.


  §.§.§    Amplitudes for B_c→ (S,A,T)P decays 

We mark  subscript   S, A, and T to denote the decay amplitudes contain a
scalar, axial-vector, and tensor charmonium in the final states, respectively.
The amplitudes from factorizable diagrams (a) and (b) in Fig. <ref> for B_c→ SP,AP decays read as

    ℱ_S   =   2√(2/3)C_Ff_Bf_Pπ M^4(r_2^2-1)∫_0^1∫_0^1dx_1dx_2 ∫_0^∞∫_0^∞b_1b_2db_1db_2ϕ_B_c(x_1)
    
          [r_2ψ_S^s(x_2,b_2)(r_b-2x_2)+ψ_S^v(x_2,b_2)(x_2-2r_b)]E_ab(t_a)h(α_e,β_a,b_1,b_2)S_t(x_2)-
    
          [ψ_S^v(x_2,b_2)(r_c+r_2^2(x_1-1))-2r_2ψ_S^s(x_2,b_2)(r_c+x_1-1)]E_ab(t_b)h(α_e,β_b,b_2,b_1)S_t(x_1),


    ℱ_A   =   -2√(2/3)C_Ff_Bf_Pπ M^4(r_2^2-1)∫_0^1∫_0^1dx_1dx_2 ∫_0^∞∫_0^∞b_1b_2db_1db_2ϕ_B_c(x_1)
    
          [r_2ψ_A^t(x_2,b_2)(r_b-2x_2)+ψ_A^L(x_2,b_2)(x_2-2r_b)]E_ab(t_a)h(α_e,β_a,b_1,b_2)S_t(x_2)+
    
          ψ_A^L(x_2,b_2)[r_c+r_2^2(x_1-1)]E_ab(t_b)h(α_e,β_b,b_2,b_1)S_t(x_1),

respectively.
The corresponding formula for nonfactorizable   diagrams  (c) and  (d) are

    ℳ_S   =   8/3C_Ff_Bπ M^4(r_2^2-1)∫_0^1∫_0^1∫_0^1dx_1dx_2dx_3
     ∫_0^∞∫_0^∞b_1b_3db_1db_3ϕ_B_c(x_1)ϕ^A_P(x_3)
    
          [ψ_S^v(x_2,b_1)(r_2^2(x_1+2x_2+x_3-2)+x_1-x_3)-r_2ψ_S^s(x_2,b_1)(x_1+x_2-1)]E_cd(t_c)h(β_c,α_e,b_3,b_1)-
    
          [ψ_S^v(x_2,b_1)(r_2^2(x_2-x_3)+2x_1+x_2+x_3-2)-r_2ψ_S^s(x_2,b_1)(x_1+x_2-1)]E_cd(t_d)h(β_d,α_e,b_3,b_1),


    ℳ_A   =   8/3C_Ff_Bπ M^4(r_2^2-1)∫_0^1∫_0^1∫_0^1dx_1dx_2dx_3
     ∫_0^∞∫_0^∞b_1b_3db_1db_3ϕ_B_c(x_1)ϕ^A_P(x_3)
    
          [ψ_A^L(x_2,b_1)(r_2^2-1)(x_1-x_3)-r_2ψ_A^t(x_2,b_1)(x_1+x_2-1)]E_cd(t_c)h(β_c,α_e,b_3,b_1)+
    
          [ψ_A^L(x_2,b_1)(r_2^2(x_2-x_3)+2x_1+x_2+x_3-2)-r_2ψ_A^t(x_2,b_1)(x_1+x_2-1)]E_cd(t_d)h(β_d,α_e,b_3,b_1),

with r_b,c=m_b,c/M. C_F=4/3 is a color factor. f_P is the decay constant of the light pseudoscalar meson, emitted from the weak vertex.
 The functions h, E and the factorization scales t_a,b,c,d can be found in <cit.>.
 The leading twist DAs of the pseudoscalar meson ϕ_P^A and
 the jet function S_t(x) come from  <cit.>.
α_e and β_a,b,c,d are the virtuality of the internal gluon and quarks, respectively. Their expressions are

    α_e   =   -[x_1-(1-x_2)r_2^2](x_1+x_2-1)M^2,
    β_a   =   [r^2_b-x_2(1+r_2^2(x_2-1))]M^2,
    β_b   =   [r_c^2+(x_1-1)(r_2^2-x_1)]M^2,
    β_c   =   -(x_1+x_2-1)[x_1-x_3+r_2^2(x_2+x_3-1)]M^2,
    β_d   =   -(x_1+x_2-1)[x_1+x_3-1+r_2^2(x_2-x_3)]M^2.

It should be stressed that    the nonlocal matrix element for the axial-vector and  scalar  charmonium meson in Eq. (<ref>)
can be  related to the vector and pseudoscalar ones <cit.>, respectively,
 by multiplying by the structure -(i)γ_5 from the left hand.
The factorization formulas (ℱ/ℳ) here and below are similar to the corresponding ones in <cit.> with some terms flipping signs.
As mentioned before, the nonlocal matrix element associating with the tensor charmonium in Eq. (<ref>)
is also analogous to the vector case,  except that the polarization vector is replaced by ϵ_∙.
Therefore  the correspondence between a tensor meson and a axial-vector
meson allows us to get the factorization formulas of B_c→ TP as


    ℱ_T=√(2/3)ℱ_A|_ψ_A^L→ψ_T, ψ_A^t→ψ_T^t,r_c→ -r_c ,
      ℳ_T=√(2/3)ℳ_A|_ψ_A^L→ψ_T, ψ_A^t→ψ_T^t.


With the functions obtained in the above, the total decay amplitudes  for the B_c→ (S,A,T)P are given by

    𝒜(B_c→ (S,A,T)P)=V^*_cbV_uq[(C_2+1/3C_1)ℱ_S,A,T+C_1ℳ_S,A,T].



  §.§.§   Amplitudes for B_c→ (S,A,T)V decays 

For B_c→ SV decays, the  decay amplitudes of factorization emission diagrams and nonfactorization emission diagrams are given as

    ℱ_S   =   2√(2/3)C_Ff_Bf_Vπ M^4√(1-r_2^2)∫_0^1∫_0^1dx_1dx_2 ∫_0^∞∫_0^∞b_1b_2db_1db_2ϕ_B_c(x_1)
    
          [r_2ψ_S^s(x_2,b_2)(r_b-2x_2)+ψ_S^v(x_2,b_2)(x_2-2r_b)]E_ab(t_a)h(α_e,β_a,b_1,b_2)S_t(x_2)-
    
          [ψ_S^v(x_2,b_2)(r_c+r_2^2(x_1-1))-2r_2ψ_S^s(x_2,b_2)(r_c+x_1-1)]E_ab(t_b)h(α_e,β_b,b_2,b_1)S_t(x_1),


    ℳ_S   =   8/3C_Ff_Bπ M^4√(1-r_2^2)∫_0^1∫_0^1∫_0^1dx_1dx_2dx_3
     ∫_0^∞∫_0^∞b_1b_3db_1db_3ϕ_B_c(x_1)ϕ_V(x_3)
    
          [ψ_S^v(x_2,b_1)(r_2^2(x_1+2x_2+x_3-2)+x_1-x_3)-r_2ψ_S^s(x_2,b_1)(x_1+x_2-1)]E_cd(t_c)h(β_c,α_e,b_3,b_1)-
    
          [ψ_S^v(x_2,b_1)(r_2^2(x_2-x_3)+2x_1+x_2+x_3-2)-r_2ψ_S^s(x_2,b_1)(x_1+x_2-1)]E_cd(t_d)h(β_d,α_e,b_3,b_1),

where f_V and ϕ_V are the decay constants and the twist-2 distribution amplitudes of the light vector mesons, respectively.
The total decay amplitudes for B_c→ SV decays are  similar to that of B_c→ SP in  Eq. (<ref>) with the replacement
f_P→ f_V, ϕ_P^A→ϕ_V.


Like vector mesons, axial-vector mesons  also carry spin degrees of freedom. Therefore,
B_c→ AV decays contain more amplitudes associated with three
different polarizations. 
We mark  superscript L, N and T to denote the contributions from longitudinal polarization, normal polarization, and transverse
polarization, respectively.


    ℱ^L_A   =   2√(2/3)C_Ff_Bf_Vπ M^4√(1-r_2^2)∫_0^1∫_0^1dx_1dx_2 ∫_0^∞∫_0^∞b_1b_2db_1db_2ϕ_B_c(x_1)
    
          [r_2ψ_A^t(x_2,b_2)(r_b-2x_2)+ψ_A^L(x_2,b_2)(x_2-2r_b)]E_ab(t_a)h(α_e,β_a,b_1,b_2)S_t(x_2)+
    
          ψ_A^L(x_2,b_2)[r_c+r_2^2(x_1-1)]E_ab(t_b)h(α_e,β_b,b_2,b_1)S_t(x_1),


    ℱ^N_A   =   2√(2/3)C_Ff_Bf_Vr_3π M^4∫_0^1∫_0^1dx_1dx_2 ∫_0^∞∫_0^∞b_1b_2db_1db_2ϕ_B_c(x_1)
    
          [ψ_A^T(x_2,b_2)(r_2^2(r_b+2-4x_2)+r_b-2)+r_2ψ_A^V(x_2,b_2)
    (r_2^2(x_2-1)+x_2+1)-4r_b]
          E_ab(t_a)h(α_e,β_a,b_1,b_2)S_t(x_2)-
    ψ_A^V(x_2,b_2)r_2[1-2r_c-2x_1+r_2^2]E_ab(t_b)h(α_e,β_b,b_2,b_1)S_t(x_1),


    ℱ^T_A   =   2√(2/3)C_Ff_Bf_Vr_3π M^4(r_2^2-1)∫_0^1∫_0^1dx_1dx_2 ∫_0^∞∫_0^∞b_1b_2db_1db_2ϕ_B_c(x_1)
    
          [ψ_A^T(x_2,b_2)(2-r_b)+r_2ψ_A^V(x_2,b_2)
    (x_2-1)]
          E_ab(t_a)h(α_e,β_a,b_1,b_2)S_t(x_2)+
    ψ_A^V(x_2,b_2)r_2E_ab(t_b)h(α_e,β_b,b_2,b_1)S_t(x_1),


    ℳ^L_A   =   8/3C_Ff_Bπ M^4√(1-r_2^2)∫_0^1∫_0^1∫_0^1dx_1dx_2dx_3
     ∫_0^∞∫_0^∞b_1b_3db_1db_3ϕ_B_c(x_1)ϕ_V(x_3)
    
          [-ψ_A^L(x_2,b_1)(r_2^2-1)(x_1-x_3)+r_2ψ_A^t(x_2,b_1)(x_1+x_2-1)]E_cd(t_c)h(β_c,α_e,b_3,b_1)-
    
          [ψ_A^L(x_2,b_1)(r_2^2(x_2-x_3)+2x_1+x_2+x_3-2)-r_2ψ_A^t(x_2,b_1)(x_1+x_2-1)]E_cd(t_d)h(β_d,α_e,b_3,b_1),


    ℳ^N_A   =   8/3C_Ff_Br_3π M^4∫_0^1∫_0^1∫_0^1dx_1dx_2dx_3
     ∫_0^∞∫_0^∞b_1b_3db_1db_3ϕ_B_c(x_1)
    
          [ϕ_V^a(x_3)ψ_A^V(x_2,b_1)(r_2^2-1)(x_2+x_3-1)
    +ϕ_V^v(x_3)ψ_A^T(x_2,b_1)(r_2^2(x_1+2x_2+x_3-2)+x_1-x_3)]
          E_cd(t_c)h(β_c,α_e,b_3,b_1)+
    [ϕ_V^a(x_3)2(r_2^2-1)(r_2(x_2-x_3)ψ_A^V(x_2,b_1)+2(x_1+x_3-1)
    ψ_A^T(x_2,b_1))
          +ϕ_V^v(x_3)(2r_2ψ_A^V(x_2,b_1)(r_2^2(x_2-x_3)+2x_1+x_2+x_3-2)
    +
          ψ_A^T(x_2,b_1)(r_2^2(1-x_1-2x_2+x_3)+1-x_1-x_3))]E_cd(t_d)h(β_d,α_e,b_3,b_1),


    ℳ^T_A   =   -8/3C_Ff_Br_3π M^4∫_0^1∫_0^1∫_0^1dx_1dx_2dx_3
     ∫_0^∞∫_0^∞b_1b_3db_1db_3ϕ_B_c(x_1)
    
          [ϕ_V^a(x_3)ψ_A^V(x_2,b_1)2r_2(r^2_2(x_2+x_3-1)+2x_1+x_2-x_3-1)
    +ϕ_V^v(x_3)ψ_A^T(x_2,b_1)(r_2^2-1)(x_1-x_3)]
          E_cd(t_c)h(β_c,α_e,b_3,b_1)+
    [2ϕ_V^a(x_3)(r_2ψ_A^V(x_2,b_1)(r_2^2(x_2-x_3)+2x_1+x_2+x_3-2)
          +2
    ψ_A^T(x_2,b_1)(r_2^2(1-x_1-2x_2+x_3)+1-x_1-x_3))
          +ϕ_V^v(x_3)ψ_A^T(x_2,b_1)
    (r_2^2-1)(x_1+x_3-1)]E_cd(t_d)h(β_d,α_e,b_3,b_1),

where ϕ_V^a and ϕ_V^v are the two twist-3 distribution amplitudes for the transverse  polarization of light vector mesons.
For B_c→ TV decays, the  decay amplitudes can be related to the axial-vector ones by making the following replacement:

    ℱ(ℳ)^L_T   =   √(2/3)ℱ(ℳ)^L_A|_ψ_A^L→ψ_T, ψ_A^t→ψ_T^t,r_c→ -r_c ,
    ℱ(ℳ)^N,T_T   =   √(1/2)ℱ(ℳ)^N,T_A|_ψ_A^T→ψ_T^T, ψ_A^V→ψ_T^V,r_c→ -r_c ,

 where the factors √(2/3) and √(1/2) come from the equivalent polarization vector
  ϵ_∙ in Eq. (<ref>) of the tensor mesons for  the longitudinal and transverse  polarizations, respectively.
  The total decay amplitudes for B_c→ (A,T)V decays can be decomposed as

    𝒜(B_c→ (A,T)V)=𝒜^L_A,T+𝒜^N_A,Tϵ_2T·ϵ_3T
    +i𝒜^T_A,Tϵ_αβρσn^αv^βϵ_2 T^ρϵ_3 T^σ,

where the three  polarization amplitudes   have the same structure as Eq. (<ref>).



§ NUMERICAL RESULTS

To proceed the numerical analysis, it is useful to summarize all of the input quantities we have used in this work.
The central values (in GeV) of the relevant meson masses and heavy quark masses are adopted as  <cit.>

    M   =   6.275,    m_b=4.8,    m_c=1.275,    m_ρ=0.775,    m_K^*=0.892.

While the masses of the P-wave charmonium have been given  in Table <ref>.
The CKM matrix-elements are set as |V_cb|=0.0405, |V_us|=0.2248, and
|V_ud|=0.97417  <cit.>.

For the decay constants of P-wave charmonium, the detailed discussions
in the nonrelativistic QCD factorization, the light-front approach and the QCD sum rules,
 could be found in Refs. <cit.>.
Here we employ the most recent updated values (in GeV) evaluated from the QCD sum rules at the scale μ=m_c <cit.>:

    f_χ_c0   =   0.0916,    f_χ_c1=0.185,    f_χ_c2=0.177,    f_h_c=0.127, 
       f^⊥_χ_c1   =   0.0875,    f^⊥_χ_c2=0.128,     f^⊥_h_c=0.133.

For the decay constants of light mesons, we use <cit.>

    f_π   =   0.131 ,   f_K=0.160 ,   f_ρ=0.209  ,
     f_K^*   =   0.217  ,     f_ρ^⊥=0.165  ,   f_K^*^⊥=0.185   GeV.

The B_c meson decay constant and lifetime  are adopted as  f_B_c=0.489 GeV <cit.>
and τ_B_c=0.507 ps <cit.>, respectively.


The branching ratios for the considered  decays in the B_c meson rest frame can be written as

    ℬ=G^2_Fτ_B_c/32π M(1-r^2_2)|𝒜|^2,

where  the decay amplitudes 𝒜 for
each channel have been given explicitly in the previous section.
When the final states involve   axial-vector/tensor charmonium states and a vector meson, the decay
amplitude can be decomposed into three components,

    |𝒜|^2=|𝒜_0|^2+|𝒜_∥|^2+|𝒜_⊥|^2,

where 𝒜_0, 𝒜_∥,𝒜_⊥ refer to the longitudinal,
parallel, and perpendicular polarization amplitudes in the transversity
basis, respectively, which are related to 𝒜^L,N,T of Eq. (<ref>) via

    𝒜_0=𝒜^L,   𝒜_∥=√(2)𝒜^N,   𝒜_⊥=√(2)𝒜^T.








Our numerical results of branching ratios for B_c→ (S,A,T)P  and B_c→ (S,A,T)V
decays are listed in Tables <ref> and <ref>, respectively.
The first kind of uncertainties is from the shape parameter
ω in the  wave function of the B_c meson and the charm-quark mass m_c.
 In the evaluation, we vary the value of ω within a 20% range and m_c = 1.275 GeV by ± 0.025 GeV.
The second error comes from the decay constants of the P-wave charmonium meson in Eq. (<ref>), which
varies 10% for error estimates.
The last one is caused  by the hard scale t located between 0.75∼ 1.25
times the invariant masses of the internal particles, which characterizes the  size of higher-order corrections to the hard amplitudes.
 It turns out that the errors are dominant by the uncertainties from the decay
 constants of P-wave charmonium meson distribution amplitudes, which can reach 20% in magnitude.
As discussed in Ref <cit.>, by using the light-cone wave function for
the B_c  meson, the theoretical uncertainty from the 
the charm-quark mass is largely reduced. It is also found that the branching ratios  are insensitive to
the hard scale, which  means the higher-order contributions can be safely neglected.
In a recent paper <cit.>, the authors claimed that the relativistic corrections to light-cone distribution
amplitudes of S-wave heavy quarkonia  are comparable with the next-to-leading order radiative corrections.
In view of this point,
we  check  the sensitivity of our results to the   squared velocity
 v^2 of the charm quark inside the P-wave charmonium  states in Eq.(<ref>).
The variation of v^2 in the range 0.25∼ 0.35 will result in the branching ratios changing only  a few percents.
This is similar to the comment in <cit.> that


 the relativistic  corrections to the Coulomb wave functions would be less significant.

In addition, the uncertainties related to the light  mesons,
such as the  decay constants and the Gegenbauer
moments shown in <cit.>, are less than 10%.
Therefore they have been neglected in our calculations.


It can be seen that
the  former four processes ( including one π or ρ meson in the  final states)
 have a relatively large branching ratios   owing to the CKM factor enhancement, while
the branching ratios of the latter four processes ( including one K or K^* meson in the final states)
are comparatively small  due to the CKM factor suppression.
Since the two type decays  have identical topology and similar kinematic properties.
In the limit of SU(3) flavor symmetry, the relative
ratios  ℛ_K/π≡ℬ(B_c→ (S,A,T)K)/ℬ(B_c→ (S,A,T)π)
and  ℛ_K^*/ρ≡ℬ(B_c→ (S,A,T)K^*)/ℬ(B_c→ (S,A,T)ρ)
are dominated by the ratio of the relevant CKM matrix
elements |V_us|^2/|V_ud|^2∼λ^2 under the naive
factorization approximation. After including  the kaon (K^*) and pion (ρ) decay constants,
one expects ℛ_K/π∼ 0.081  and ℛ_K^*/ρ∼ 0.057.
From Tables <ref> and <ref>,
our  predictions for ℛ_K/π corresponding to various P-wave charmonium states  lie in the range   0.075 to 0.080, while ℛ_K^*/ρ  
is in the range 0.057 to 0.064,  both are very close to the above expected values.
It means that the dominant contributions to the branching ratios come from the factorizable
topology, while the nonfactorizable contribution is suppressed by the Wilson coefficient C_1
[see Eq. (<ref>)]. 


One can see some  interesting hierarchical relations among these branching ratios in our predictions. For example,
branching ratios for decays involving  pseudoscalar mesons
in the final state are smaller than their vector partners for the same flavor content.
This is partially due to the pseudoscalar meson decay constant is usual smaller than the vector ones.
Furthermore, since the B_c meson is a spinless particle,
 according to the  angular momentum conservation,
only one partial wave contribute to the B_c→ SP, AP, TP, SV decays,
while in the AV, TV modes,  three partial waves are simultaneously allowed, resulting in the larger branching ratios. 
For those channels with the  same light meson and different P-wave charmonium mesons
 in the  final states, we have the following hierarchy pattern:

    ℬ(B_c→χ_c2P)>ℬ(B_c→χ_c0P)>
    ℬ(B_c→χ_c1P)∼ℬ(B_c→ h_c P).

As discussed   in Ref <cit.>, the branching ratio of B_c→η_c(2S) π is enhanced by the
the twist-3 distribution amplitude from Fig. 3(b). Nevertheless, 
this contribution vanishes   for  B_c→ψ(2S)π decay
since the Lorentz structure of the vector charmonium wave functions is different from the pseudoscalar case.
As mentioned in the previous section,
because of the corresponding relation between a pseudoscalar (vector) and a scalar (axial-vector) charmonium,
the similar situation also exists in this work.
The twist-3 distribution amplitude  from Fig <ref>(b) also give the dominant contribution to the
B_c→χ_c0P decays, while for other channels, the dominant contribution still come from the twist-2 ones.
Because the strong interference between the twist-2 and twist-3 contributions   is constructive in B_c→χ_c0P,
we   have a large branching ratio for this mode.
One can see that the dominant twist-2 contributions for B_c→ (χ_c1,h_c) P
 are suppressed by a factor of r_c-r^2_2 given in Eq. (<ref>),
whereas this  suppression  is absent in the case of B_c→χ_c2 P
 due to the r_c term flipping sign [see Eq. (<ref>)].
 This explains why  B_c→χ_c2P has a rate greater than χ_c1P and h_cP.
Of course,  this is only a rough estimate on the magnitudes, the branching ratios also have been related
to the decay constants and the distribution amplitudes of the various P-wave charmonium mesons.
The  relations between the decay  constants f_χ_c2>f_χ_c0 in Eq. (<ref>) implies that
ℬ(B_c→χ_c2P)>ℬ(B_c→χ_c0P).
The similar pattern also occurs B_c→ (A,S,T)V decays; see Table <ref>.


As mentioned in the Introduction,  many   other work have performed  a systematic study on
the P-wave charmonium decays of B_c mesons. Various approaches such as
several relativistic and nonrelativistic quark models <cit.>,
the  sum rules of QCD <cit.>,
the improved Bethe-Salpeter approach <cit.>, the Isgur-Scora-Grinstein-Wise II model <cit.>,
 and the nonrelativistic QCD effective theory <cit.>
have been used to calculate the branching ratios.
For the sake of comparison, we  briefly list the obtained  theoretical results in Tables <ref> and <ref>.
One finds that some of  the results  given by   different models are roughly comparable.
For example, our theoretical predictions on those decays involving the h_c meson in the final state
are of the same order of magnitude as observed in <cit.>.
The branching ratios of B_c→χ_c1π and
 B_c→ h_cπ evaluated by N. Sharma et al. <cit.>
 are 7× 10^-4 and 6× 10^-4, respectively,  which also match well with our results.
 In a very recent paper <cit.>, the author predicted the branching ratio
ℬ(B_c^+→χ_c0π^+)=1.22× 10^-3, which is comparable
to our prediction. Of course, some predicted values are quite a spread in various models.
The predictions in Ref <cit.> are typically larger excepted for B_c^+→χ_c1π^+.
Previously, G. L. Castro et al. <cit.> studied the nonleptonic decays of the B_c into tensor mesons using
the factorization hypothesis. They predict ℬ(B_c→χ_c2π),  ℬ(B_c→χ_c2K),
ℬ(B_c→χ_c2ρ),  and ℬ(B_c→χ_c2K^*), as 7.5× 10^-5,
5.49× 10^-6, 2.38× 10^-4, and 1.33× 10^-5, respectively,
which are considerably smaller than our results as well as most of other model calculations.
Our results  for a final K^(*)  are also larger than those of other approaches.
The disagreement in the predictions may be attributed to the
different values of the form factors used in these approaches.
Experimental investigations on these decays may be used to test theoretical methods according to their predictions.



On the experimental side, so far only the evidence  for the decay
B^+_c→χ_c0π^+ is found at 4.0 σ  significance by the LHCb Collaboration <cit.>.
The ratio of production cross sections of the B^+_c and B^+ mesons times branching fractions is measured to be
σ_B_c^+/σ_B^+×ℬ(B^+_c→χ_c0π^+)
=(9.8^+3.4_-3.0(stat)± 0.8(syst))× 10^-6 <cit.>.
As a cross-check,  the   cross section  ratio σ_B_c^+/σ_B^+ can be extracted from 
another charmonium mode,  σ_B_c^+/σ_B^+×ℬ(B^+_c→ J/ψπ^+)/ℬ(B^+→ J/ψ K^+)=(0.683±0.018±0.009)%  measured by the
LHCb Collaboration <cit.>. The branching ratio ℬ(B^+→ J/ψ K^+), determined from the world average value,
is (1.026±0.031)×10^-3 <cit.>.
If we use our previous PQCD calculation  ℬ(B^+_c→ J/ψπ^+)=(2.33^+0.81_-0.61)× 10^-3 <cit.>,
where all errors are combined in quadrature, as an input, the ratio  σ_B_c^+/σ_B^+ is in the region of
(2.2∼ 4.1)× 10^-3. 
Combined with the prediction on  ℬ(B^+_c→χ_c0π^+) in Table <ref>, we obtain the range  σ_B_c^+/σ_B^+×ℬ(B^+_c→χ_c0π^+)=(2.6∼ 8.2)× 10^-6,
which is consistent with the LHCb data  with one sigma errors.






Turning to the polarizations 
for B_c → AV, TV decays.
 We usually define five observables corresponding to three polarization fractions f_λ(λ=0,∥,⊥),
and two relative phases ϕ_∥, ϕ_⊥, where
 
    f_λ=|𝒜_λ|^2/|𝒜_0|^2+|𝒜_∥|^2+|𝒜_⊥|^2,  ϕ_∥,⊥=arg𝒜_∥,⊥/𝒜_0,

with normalization such that ∑_λf_λ=1.
The results for the   polarization fractions
and their relative phases
are displayed in Table <ref>, where the sources of the errors in the
 numerical estimates have the same origin as in the discussion of the branching ratios in Table <ref>.
 It can be observed that both the polarization fractions and the phases 
 are relatively stable with respect to the variations of hadronic parameters, the decay constants and the hard scale,
 and therefore they serve as  good quantities to test the standard model.
Several remarks are given in order. First, the contributions to the  branching
ratios mainly arise from the longitudinal polarizations 
because of the relation f_0≫ f_∥∼ f_⊥, which is expected
from the power counting rules.
For example,
the longitudinal parts of B_c→ TV  decays occupy over 90%,
which are very similar to the case of B_c→ J/ψ V <cit.>.
However,  the longitudinal polarizations of B_c→χ_c1V are relative smaller (∼60%)
compared to that of B_c→ h_c V.
As mentioned before, owing to the G-parity, the  distribution amplitudes  for χ_c1 and
h_c  mesons exhibit the different asymptotic behaviors (see Eqs. (<ref>) and (<ref>)).
If we use the h_c distribution amplitudes for calculation,  the resultant predictions  f_0(χ_c1V)
 can be increase to around 90%. 
 Besides, the longitudinal and transverse decay constants in the two axial-vector mesons can also contribute to
 different polarizations.
Second, for B_c→ (A,T)ρ and B_c→ (A,T)K^* decays,
both   have similar magnitudes and phases of the amplitudes, which
suggests the SU(3) breaking effect between them is small.
Last, the predicted relative phases deviations from π
indicate the existence of the still unknown final-state interaction.
 However, the magnitudes and phases of the two transverse
amplitudes 𝒜_∥ and 𝒜_⊥ are roughly equal,
which is expected from analyses based on quark-helicity conservation <cit.>.

These results and findings will be further tested by the LHCb and
Belle-II experiments in the near future.





§  CONCLUSION


The two-body  B_c meson decays to a P-wave charmonium state (χ_c0,χ_c1,χ_c2,h_c) and a light (π,K,ρ,K^*) meson
 are systematically  analysed within the  perturbative QCD approach.
Our predictions for the branching ratios are summarized in Tables <ref> and <ref> and compared with other theoretical
results. Overall, the  predicted branching ratios from different theoretical models  have a relative big spread.
The upcoming experimental measurements of the corresponding decay rates can examine various theoretical approaches.
Based on our estimations, the dominating decay mode of the concerned processes is B_c→χ_c2ρ
with predicted  branching ratios of 1.6%, which should   be
accessible experimentally at high-luminosity hadron colliders.
We also estimate the polarization contributions in B_c→ (χ_c1,c2,h_c)V decays.
As expected, based on the factorization assumption,
the longitudinal polarization dominates and
the transverse polarizations are of the same size.


We also discussed theoretical uncertainties arising from
the hadronic parameters in B_c meson wave function, the decay constants of  charmonium  states and the hard scale t.
The branching ratios suffer a large error from the decay constants, whereas the polarization observables are less sensitive
to these parameters. The obtained results can be confronted to the experimental data in the future.


I would like to acknowledge Hsiang-nan Li and Yu-Ming
Wang for helpful discussions. This work is supported in
part by the National Natural Science Foundation of China
under Grants No. 11547020 and No. 11605060 and, in part,
by the Program for the Top Young Innovative Talents of
Higher Learning Institutions of Hebei Educational
Committee under Grant No. BJ2016041.








§ DETAILS FOR DERIVING THE P-WAVE CHARMONIUM DAS

Starting with the  momentum-space radial wave function which can be written as  the Fourier transform of the position-space expression ψ_nlm(r⃗)

    ψ(k)=∫_-∞^∞ψ_nlm(r⃗)e^-ik⃗·r⃗dr⃗,

where n, l, and m stand for main, orbital, and magnetic quantum  numbers, respectively.
In above equation, the first term  ψ_nlm(r⃗) is known to be separated into R_nl(r)Y_lm(θ,φ)
 in the spherical coordinates (r, θ, φ),
while  the second exponential term  in the plane wave expansion can be written as

    e^-ik⃗·r⃗=e^-ik r cosθ=∑_l'=0^∞√(4π(2l'+1))(-i)^l'j_l'(kr)Y_l'0(θ,0),

with j_l'(kr)  the spherical Bessel function.
We then write  Eq.(<ref>)  as

    ψ(k)= √(4π(2l+1))(-i)^l∫_0^∞j_l(kr) R_nl(r) r^2 dr,

where  the orthogonality property 
∫_0^π∫_0^2πY_lmY_l'0sinθ dθ dφ=δ_ll'δ_m0
have been used.

For the P-wave states n=2 and l=1, employing the spherical Bessel function j_1(kr)=sin(kr)-krcos(kr)/(kr)^2
and the radial wave function  for a Coulomb Potential R_21(r)∝ r e^-q_B r/2,
 the integral of Eq. (<ref>) evaluates to

    ψ(k) ∝k q_B/(4k^2+q_B^2)^3,

where q_B is the Bohr momentum.
Note that the above expression  is in contrast to   Eq.(47) in <cit.>.
We argue that  the spherical harmonics function for P-wave states is dependent on the angle θ, which should  contribute to the
integral in Eq. (<ref>).
In particular,  Eq.(<ref>) is almost the same as M. Beneke's calculation in Ref. <cit.> (see Eq.(45)),
except for a constant term which can be absorbed in the redefinition of the wave function of the P-wave charmonium.

Following much the same procedure as described in Refs <cit.>,
we obtain the heavy quarkonium DA which is dependent on the charm quark momentum fraction x after integrating the transverse
momentum k_T,

    Φ(x)∼∫ d^2k_T ψ(x,k_T)∝ x(1-x){√(x(1-x)(1-4x(1-x))^3)/[1-4x(1-x)(1-v^2/4)]^2},

where v=q_B/m_c is the charm quark velocity. In the numerical calculation, we take v^2=0.3 and neglect the v^2 term
 in the numerator  <cit.>.
 As mentioned in Eq. (<ref>), we propose the P-wave charmonium states DAs  as
ψ(x)∝Φ_asy(x)𝒯(x)
with
 
    𝒯(x)={√(x(1-x)(1-4x(1-x))^3)/[1-4x(1-x)(1-v^2/4)]^2}^1-v^2,

where the power 1-v^2 denotes the small relativistic corrections to the Coulomb wave functions <cit.>.




99
pdg2016
C. Patrignani et al. (Particle Data Group), Chin. Phys. C 40, 100001 (2016).
prd74051103
R. Kumar   et al. (Belle Collaboration), Phys.Rev. D 74, 051103(R) (2006).
prd78091104
R. Kumar et al. (Belle Collaboration), Phys.Rev. D  78, 091104 (2008).
prl88031802
K. Abe et al. (Belle Collaboration), Phys. Rev. Lett. 88, 031802 (2002).
prd78091101
B. Aubert et al. (BABAR Collaboration), Phys. Rev. D 78,  091101(R) (2008).
prl89011803
K. Abe et al. (Belle Collaboration), Phys. Rev. Lett. 89, 011803 (2002).
prl94141801
B. Aubert et al. (BABAR Collaboration), Phys. Rev. Lett. 94, 141801 (2005).
prl102132001
B. Aubert et al. (BABAR Collaboration), Phys. Rev. Lett. 102, 132001 (2009).
plb643155
N. Soni  et al. (Belle Collaboration), Phys. Lett. B 643, 155 (2006).
prl107091803
V. Bhardwaj et al. (Belle Collaboration), Phys. Rev. Lett. 107,  091803 (2011).

prd74012007
F. Fang, et al. (Belle Collaboration), Phys. Rev. D   74, 012007 (2006).
prd78012006
B. Aubert et al. (BABAR Collaboration), Phys. Rev. D   78, 012006 (2008).
plb874663
R. Aaij et al. (LHCb Collaboration),  Nucl. Phys.  B874, 663 (2013).


prd85052003
J. P. Lees et al. (BABAR Collaboration), Phys. Rev. D 85, 052003 (2012).
prd93052016
V. Bhardwaj et al. (Belle Collaboration), Phys. Rev. D  93, 052016 (2016).
prl119062001
R. Aaij et al. (LHCb Collaboration),  Phys. Rev. Lett. 119, 062001 (2017).
prd94091102
R. Aaij et al. (LHCb Collaboration),  Phys. Rev. D 94, 091102(R) (2016).



prd493399
C. H. Chang and Y. Q. Chen, Phys. Rev. D 49, 3399 (1994).
prd564133
J.F. Liu, K.T. Chao, Phys. Rev. D 56, 4133 (1997).
prd62014019
 A. Abd El-Hady, J. H. Munoz, and J. P. Vary, Phys. Rev. D 62, 014019 (2000).
prd61034012
P. Colangelo and F. De Fazio, Phys. Rev. D 61, 034012 (2000).
npb585353
V. V. Kiselev, A. E. Kovalsky, and A. K. Likhoded, Nucl. Phys. B585, 353 (2000).
prd68094020
D. Ebert, R. N. Faustov, and V. O. Galkin, Phys. Rev. D 68, 094020 (2003).


jhep06015
H. Fu, Y. Jiang, C. S. Kim, and G. L. Wang, J. High Energy Phys. 06, 015 (2011).
prd86094028

 Sk. Naimuddin, S. Kar, M. Priyadarsini, N. Barik, and P. C. Dash, Phys. Rev. D 86, 094028 (2012).
prd88094014
S. Kar, P. C. Dash, M. Priyadarsini,  Sk. Naimuddin, and N. Barik, Phys. Rev. D 88, 094014 (2013).
prd89034008
C.-F. Qiao, P. Sun, D. Yang, and R.-L. Zhu, Phys. Rev. D 89, 034008 (2014).
prd89017501
H.W. Ke, T. Liu, and X.Q. Li, Phys. Rev. D 89, 017501 (2014).
scp58071001
C.H. Chang, H.F. Fu, G.L. Wang, and J.M. Zhang,  Sci China-Phys Mech Astron 58, 071001 (2015).




prd82034019
D. Ebert, R. N. Faustov, and V. O. Galkin, Phys. Rev. D  82, 034019 (2010).

prd74074008
E. Hernández, J. Nieves, and J. M. Verde-Velasco, Phys. Rev. D  74, 074008 (2006).
prd71094006
M. A. Ivanov, J.G. Körner, and P. Santorelli, Phys. Rev. D 71, 094006 (2005).
prd73054024
M. A. Ivanov, J.G. Körner, and P. Santorelli, Phys. Rev. D 73, 054024 (2006).
prd77054003
Y. M. Wang, C. D. Lü, Phys. Rev.D 77, 054003 (2008).
prd65014017
C. H. Chang, Y. Q. Chen, G. L. Wang, and H. S. Zong,
Phys. Rev. D 65, 014017 (2001); Commun. Theor. Phys.
35, 395 (2001).
jpg28595
V.V. Kiselev, O. N. Pakhomova, and V. A. Saleev, J. Phys. G 28, 595 (2002).
jpg39015009
Zhi-hui Wang, Guo-Li Wang, and Chao-Hsi Chang, J. Phys. G 39, 015009 (2012).
prd87034004
R. Dhir and C. S. Kim, Phys. Rev. D 87, 034004 (2013).
prd79114018
X.-X. Wang, W. Wang, and C. D. Lü, Phys. Rev. D 79, 114018 (2009).

epjc76454
Q. Li, T.-h. Wang, Y. Jiang, H. Yuan, and  Guo-Li Wang, Eur. Phys. J. C 76, 454 (2016).
171007011
R. Zhu, arXiv: arXiv:1710.07011.
rui1
Z. Rui, Z.T. Zou, Phys. Rev. D 90, 114030 (2014).
rui2
Z. Rui, W.F. Wang, G.X. Wang, L.H. Song, C.D. Lü, Eur. Phys. J. C 75, 293 (2015).


plb612215
A.E. Bondar, V.L. Chernyak, Phys. Lett. B 612, 215 (2005).

rmp681125
G. Buchalla, A. J. Buras and M. E. Lautenbacher, Rev. Mod. Phys. 68, 1125 (1996).
prd82054019
H.Y. Cheng, Y. Koike, and K. C. Yang, Phys. Rev. D 82, 054019 (2010).
prd83034001
H.Y. Cheng and K.-C. Yang, Phys. Rev. D 83, 034001 (2011).
npb591313
M. Beneke, G. Buchalla, M. Neubert, C.T. Sachrajda, Nucl. Phys. B591, 313 (2000).
prd55272
 A.G. Grozin, M. Neubert, Phys. Rev. D 55, 272 (1997).
 npb5923
 M. Beneke, T. Feldmann, Nucl. Phys. B592, 3 (2001).
 epjc28515
 C.D. Lü, M.Z. Yang, Eur. Phys. J. C 28, 515 (2003).
 prd76074018
 A. Ali, G. Kramer,Y. Li, C.D. Lü, Y.L. Shen, W.Wang, Y.M.Wang,
Phys. Rev. D 76, 074018 (2007).
 prd89114019
J. Sun, Y. Yang, Q. Chang, G. Lu, Phys. Rev. D 89, 114019 (2014).
 prd71114008
C.-H. Chen, H.-N. Li, Phys. Rev. D 71, 114008 (2005).

npb776187
Kwei-Chou Yang, Nucl. Phys. B776, 187 (2007).
jhep10108
K. C. Yang, J. High Energy Phys. 10, 108 (2005).
jhep06121
X. P. Wang and D. Yang, J. High Energy Phys. 06, 121 (2014).


pr112173
V.L. Chernyak and A.R. Zhitnitsky, Phys. Rep. 112, 173 (1984).
jhep06067
M. Diehl and G. Hiller, J. High Energy Phys. 06, 067 (2001).
prd67094011
Chuan-Hung Chen, Phys. Rev. D 67, 094011 (2003).
prd74114010
Wei Wang, Yue-Long Shen, Ying Li, and Cai-Dian Lü, Phys. Rev. D 74, 114010 (2006).
prd73014017
H.Y. Cheng,  Chun-Khiang Chua,  and K. C. Yang, Phys. Rev. D 73, 014017 (2006).

prd83014008
Wei Wang, Phys Rev D 83, 014008 (2011).
ball
V. M. Braun and I. E. Filyanov, Z. Phys. C 48, 239 (1990);
P. Ball, V. M. Braun, Y. Koike, and K. Tanaka, Nucl. Phys. B529, 323 (1998);
P. Ball, J. High Energy Phys. 01, 010 (1999).
jhep10074
C. W. Hwang, J. High Energy Phys. 10, 074 (2009).
prd79074004
V.V. Braguta, A. K. Likhoded, and A.V. Luchinsky, Phys. Rev. D 79, 074004 (2009).
prd96014026
M. A. Olpak, A. Ozpineci, and V. Tanriverdi, Phys. Rev. D 96, 014026 (2017).
epjc76564
Z. Rui, H. Li, G.X.Wang, Y. Xiao, Eur. Phys. J. C 76, 564 (2016).
jhep12012
Wei Wang, Ji Xu, D. Yang, S. Zhao, J. High Energy Phys. 12, 012 (2017).




jpg37075013
N. Sharma, R. Dhir and R C Verma, J. Phys. G  37, 075013 (2010).

170900550
B. Mohammadi, J. High Energy Phys. 12, 019 (2017).
jpg282241
G. L. Castro, H. B. Mayorga, and J. H. Muñoz, J. Phys. G 28, 2241 (2002).
prl114132001
R. Aaij et al. (LHCb Collaboration), Phys. Rev. Lett. 114, 132001 (2015).
zpc1269
A. Ali, J. G. Körner, G. Kramer, and J. Willrodt, Z. Phys. C 1, 269 (1979).
prd64117503
M. Suzuki, Phys. Rev. D 64, 117503 (2001).

npb811155
M. Beneke, and L. Vernazza,  Nucl. Phys.  B811, 155 (2009).






