$B^{\ast}$ ${\to}$ $\bar{D}D$ decays with perturbative QCD approach
Junfeng Sun Haiyan Li Yueling Yang Na Wang Qin Chang Gongru Lu

  
  
  Institute of Particle and Nuclear Physics,
              Henan Normal University, Xinxiang 453007, China
  
  Institute of Particle and Nuclear Physics,
              Henan Normal University, Xinxiang 453007, China
  
  Institute of Particle and Nuclear Physics,
              Henan Normal University, Xinxiang 453007, China
  
  Institute of Particle and Nuclear Physics,
              Henan Normal University, Xinxiang 453007, China
  
  Institute of Particle and Nuclear Physics,
              Henan Normal University, Xinxiang 453007, China
  
  Institute of Particle and Nuclear Physics,
              Henan Normal University, Xinxiang 453007, China
  
  
  The nonleptonic two-body B^∗ → DD
  weak decays are studied phenomenologically with the perturbative
  QCD factorization approach.
  It is found that the B_s^∗0 → D_s^-D_s^+,
  B_d^∗0 → D_d^-D_s^+, and
  B_u^∗+ → D_u^0D_s^+
  decays have branching ratios ≳ 10^-9, and
  might be promisingly measurable at the running LHC and
  forthcoming SuperKEKB experiments in the future.
  
  12.15.Ji 12.39.St 13.25.Hw 14.40.Nd
  B^∗ → DD decays with perturbative QCD approach
    Gongru Lu
    Received: date / Accepted: date
==============================================



  
  

§ INTRODUCTION

  
  The B_q^∗ mesons, consisting of b̅q pair with
  q = u, d and s, are spin-triplet ground vector states
  with definite spin-parity quantum numbers of J^P =
  1^- <cit.>.
  Because the mass splittings m_B_q^∗ - m_B_q
  ≲ 50 MeV <cit.> are much smaller than the mass of the
  lightest pion meson, the B_q^∗ meson decays
  dominantly into the ground pseudoscalar B_q meson through
  the electromagnetic interaction.
  Besides, the B_q^∗ mesons can also decay via the
  bottom-changing transition induced by the weak interaction
  within the standard model (SM).
  Because of the strong phase-space suppression from their
  dominant magnetic dipole (M1) transition B_q^∗
  → B_qγ,
  the lifetime of the B_q^∗ meson is of the order of
  10^-17 second or less,
  which, in general, is too short to enable the B_q^∗
  meson to experience the weak disintegration <cit.>.
  The B_q^∗ weak decays have not actually attracted
  much attention yet.
  Until now, there has been no experimental report and few
  theoretical works concentrating on the B_q^∗ weak
  decay, subject to the relatively inadequate statistics on
  the B_q^∗ mesons.
  Fortunately, the high luminosities and large production
  rates at LHC and the forthcoming SuperKEKB are
  promising, and the rapid accumulation of more and more
  B_q^∗ data samples is expected to be possible.
  Some B_q^∗ weak decay modes might be detected
  and investigated in the future, which undoubtedly makes the
  B_q^∗ mesons another a vibrant arena for testing
  the Cabibbo-Kabayashi-Maskawa (CKM) picture for CP-violating
  phenomena, examining our comprehension of the underlying
  dynamical factorization mechanism, and so on.
  In addition, heavy quark symmetry relates hadronic transition
  matrix elements (HTME) of the B_q^∗ and B_q weak decays.
  The interplay between the B_q^∗ and B_q weak decays
  could prove useful information to overconstraint parameters in
  the SM, and might shed some fresh light on various anomalies in
  B decays.

  The purely leptonic decays B_q^∗ → ℓ^+ℓ^-
  induced by the flavor-changing neutral currents have been
  studied recently in the SM <cit.>.
  The semileptonic and nonleptonic B_q^∗ decays have
  been investigated also in the SM <cit.>,
  where the transition form factors are evaluated with the
  Wirbel-Stech-Bauer approach <cit.>, and the nonfactorizable
  corrections to HTME are considered <cit.>
  based on the collinear-based and QCD-improved
  factorization (QCDF) approach <cit.>.
  In this paper, we will study the nonleptonic B_q^∗
  decay into the pseudoscalar charmed-meson pair DD with the
  perturbative QCD factorization (pQCD) approach <cit.>,
  just to provide a ready reference for the future experimental
  research.
  In addition, as is well known, the production ratio for
  the B_q^∗ meson is comparable with that for the B_q
  meson (see Table <ref>), the B_q^∗ and
  B_q mesons have nearly equal mass.
  Hence, the study of the B_q^∗ → DD
  decays will undoubtedly be helpful to the experimental background
  analysis on the B_q → DD decays.

  This paper is organized as follows.
  The theoretical framework and amplitudes for B_q^∗
  → DD decays with pQCD approach are given
  in section <ref>.
  Section <ref> is devoted to numerical results and discussion.
  The final section is a summary.

  
  

§ THEORETICAL FRAMEWORK

  
  
  

 §.§ The effective Hamiltonian

  
  The effective Hamiltonian describing the B_q^∗ →
  DD weak decay is written as <cit.>
 
   
    H_ eff  = G_F/√(2)∑_q=d,s{∑_p=u,c V_pb^∗ V_pq∑_i=1^2 C_i(μ) Q_i(μ)
       -V_tb^∗ V_tq ∑_k=3^10 C_k(μ) Q_k(μ) }
       +  h.c.,

 
  where the Fermi coupling constant G_F ≃
  1.166×10^-5  GeV^-2 <cit.>;
  V_pb^∗V_pq and V_tb^∗V_tq are the CKM factors;
  The scale μ factorizes the physical contributions into two parts:
  the Wilson coefficients C_i and the local four-fermion operators Q_i.
  The operators are defined as follows.
 
    
    Q_1   =   
       ( b̅_α  p_α )_V-A 
       ( p̅_β  q_β )_V-A,
                    
        Q_2 =   ( b̅_α  p_β )_V-A 
       ( p̅_β  q_α )_V-A, 
    
        Q_3   =   ∑_q^' 
       ( b̅_α  q_α )_V-A 
       ( q̅^'_β  q^'_β )_V-A,
                 
       Q_4 = ∑_q^' 
       ( b̅_α  q_β )_V-A 
       ( q̅^'_β  q^'_α )_V-A, 
    
        Q_5   =   ∑_q^' 
       ( b̅_α  q_α )_V-A 
       ( q̅^'_β  q^'_β )_V+A,
                 
        Q_6 = ∑_q^' 
       ( b̅_α  q_β )_V-A 
       ( q̅^'_β  q^'_α )_V+A, 
    
        Q_7   =   ∑_q^' 3/2 Q_q^' 
       ( b̅_α  q_α )_V-A 
       ( q̅^'_β  q^'_β )_V+A,
               Q_8 = ∑_q^' 3/2 Q_q^' 
       ( b̅_α  q_β )_V-A 
       ( q̅^'_β  q^'_α )_V+A, 
    
        Q_9   =   ∑_q^' 3/2 Q_q^' 
       ( b̅_α  q_α )_V-A 
       ( q̅^'_β  q^'_β )_V-A,
           
       Q_10 = ∑_q^' 3/2 Q_q^' 
       ( b̅_α  q_β )_V-A 
       ( q̅^'_β  q^'_α )_V-A,

 
  where Q_1,2 are tree operators arising from the W-boson
  exchange; Q_3,⋯,6 and Q_7,⋯,10 are called
  the QCD and electroweak penguin operators, respectively;
  (q̅_1 q_2)_V±A ≡
  q̅_1 γ_μ(1±γ_5) q_2;
  α and β are color indices;
  q^' denotes all the active quarks at the scale of
  O(m_b), i.e., q^' = u, d, c, s, b;
  and Q_q^' is the electric charge of quark q^'
  in the unit of |e|.

  The Wilson coefficients C_i(μ), which summarize the
  physical contributions above the scale of μ, have
  been properly calculated at the next-to-leading order with
  the renormalization group equation assisted perturbation
  theory <cit.>.
  Due to the presence of long-distance QCD effects and the
  entanglement of nonperturbative and perturbative shares,
  the main obstacle to evaluate the B_q^(∗) weak decays
  is the treatment of physical contributions below the
  scale of μ which are included in the HTME of local operators.

  
  

 §.§ Hadronic matrix elements

  
  Some phenomenological models have recently been developed to improve
  the sketchy treatment with naive factorization scheme <cit.>.
  These models are generally based on the Lepage-Brodsky approach
  <cit.> and some power counting rules in parameters of
  α_s and Λ_ QCD/m_Q (where α_s
  is the strong coupling, Λ_ QCD is the QCD characteristic
  scale, and m_Q is the mass of a heavy quark), and express the
  HTME as a convolution integral of universal wave functions
  and hard scattering subamplitudes, such as the QCDF approach
  <cit.>, pQCD approach <cit.>,
  the soft and collinear effective theory <cit.>,
  and so on, which have been extensively employed in the interpretation
  of the B weak decays.
  To wipe out the endpoint singularities appearing in the collinear
  approximation <cit.>,
  it is suggested by the pQCD approach <cit.>
  that the transverse momentum k_T of valence quarks should be
  retaken, and a Sudakov factor should be introduced for each wave
  function to further suppress the soft contributions and make the
  hard scattering more perturbative.
  Finally, a decay amplitude is written as a multidimensional
  integral of many parts <cit.>, including the Wilson
  coefficients C_i, the heavy quark decay subamplitudes
  H, and the universal wave functions Φ,
  
  
    A ∼ ∑_i∫∏_jdk_j 
      C_i(t)  H_i(t,k_j) Φ_j(k_j) e^-S_j,

  
  where t is a typical scale; k_j is the momentum of a valence
  quark; e^-S_j is a Sudakov factor.

  
  

 §.§ Kinematic variables

  
  The light-cone variables in the rest frame of the B^∗ meson
  are defined as follows.
  
  
    p_B^∗  =  p_1  = m_1/√(2)(1,1,0)
      ,

  
  
    p_D  =  p_2  =  (p_2^+,p_2^-,0)
      ,

  
  
    p_D  =  p_3  =  (p_3^-,p_3^+,0)
      ,

  
  
    k_i  =  x_i p_i+(0,0,k⃗_iT)
      ,

  
  
    p_i^±  =  (E_i ± p)/√(2),

  
  
    ϵ_B^∗^∥  = 1/√(2)(-1,1,0)
      ,

  
  
    s  =  2 p_2·p_3,

  
  
    t  =  2 p_1·p_2  = 2 m_1 E_2,

  
  
    u  =  2 p_1·p_3  = 2 m_1 E_3,

  
  
    s t +s u-t u -4 m_1^2 p^2 = 0
      ,

  
  where the subscripts i = 1, 2, 3 of variables (energy E_i,
  momentum p_i and mass m_i) correspond to B^∗, D
  and D mesons, respectively; k_i is the momentum of the valence quark
  with the longitudinal momentum fraction x_i and the transverse momentum
  k⃗_iT; ϵ_B^∗^∥ is the longitudinal
  polarization vector; p is the momentum of the final
  states; s, t and u are the Lorentz invariant parameters.
  The notation is displayed in Fig.<ref>.

  
  

 §.§ Wave functions

  
  Wave functions are the basic input parameters with the pQCD approach.
  Although wave functions contain soft and nonperturbative contributions,
  they are universal, i.e., process independent.
  Wave functions and/or distribution amplitudes (DAs) determined
  by nonperturbative methods or extracted from data, can be employed
  here to make predictions.

  Following the notations in Refs. <cit.>,
  HTME of the diquark operators is defined as
  
  
    ⟨0|b̅_i(0)q_j(z)|B^∗(p,ϵ^∥)⟩ 
     = f_B^∗/4∫d^4k e^-ik·z{ϵ̸^∥ [
      m_B^∗ Φ_B^∗^v(k) -
      p̸ Φ_B^∗^t(k) ] }_ji,

  
  
  
    ⟨D(p)|c_i(0)q̅_j(z)|0⟩ 
     = i f_D/4∫d^4k e^+ik·z {γ_5 [ p̸ Φ_D̅^a(k)
      +m_D̅ Φ_D̅^p(k) ] }_ji,

  
  
  where f_B^∗ and f_D are decay constants; the
  wave functions Φ_B^∗^v and Φ_D̅^a
  are twist-2; and Φ_B^∗^t and Φ_D̅^p
  are twist-3.
  Due to the kinematic relation ϵ_B^∗^⊥·p_i
  = 0, the transversely polarized B^∗ meson contributes
  nothing to the amplitudes for the B^∗ → DD decays.
  The expressions for DAs of the B^∗ and D mesons are
  <cit.>
  
   
    ϕ_B^∗^v(x) = A  x x̅ exp{ -x̅ m_q^2+x m_b^2/8 ω_1^2 x x̅},

  
   
    ϕ_B^∗^t(x) = B  (x̅-x)^2 exp{ -x̅ m_q^2+x m_b^2/8 ω_1^2 x x̅},

  
  
   
    ϕ_D̅^a(x) = C  x x̅ exp{ -x̅ m_q^2+x m_c^2/8 ω_2^2 x x̅},

  
   
    ϕ_D̅^p(x) = D exp{ -x̅ m_q^2+x m_c^2/8 ω_2^2 x x̅},

  
   where x and x̅ (≡ 1 - x) are the
   parton momentum fractions;
   ω_i determines the average transverse momentum
   of partons and ω_i ≃ m_i α_s(m_i);
   parameters A, B, C, D are the normalization coefficients
   satisfying the conditions
  
   
    ∫_0^1dx ϕ_B^∗^v,t(x)=1
       ,

  
   
    ∫_0^1dx ϕ_D̅^a,p(x) =1
       .

  

  In fact, there are many phenomenological models of DAs for the
  charmed meson, for example, some of them have been listed by
  Eq.(30) in  Ref.<cit.>.
  One of the favorable models from the experimental data within the pQCD
  framework has the expression <cit.>
  
   
    ϕ_D(x) = 6 xx̅ {1+C_D(1-2x) },

  
  where parameter C_D = 0.5 for the D_u,d meson,
  and C_D = 0.4 for the D_s meson.
  In the actual calculation <cit.>,
  there is no distinction between twist-2 and twist-3 DAs.

  The shape lines of the normalized DAs ϕ_B^∗^v,t(x)
  and ϕ_D^a,p(x) are illustrated in Fig.<ref>.
  It is clearly seen from Fig.<ref> that
  (1) the shape lines of DAs in Eqs.(<ref>)-(<ref>) have a
  broad peak in the small x regions, which is generally consistent
  with an ansatz in which a light quark carries fewer parton momentum
  fractions than a heavy quark in a heavy-light system.
  (2) Due to the suppression from exponential functions, the DAs of
  Eqs.(<ref>)-(<ref>) converge quickly to zero
  at endpoint x, x̅ → 0, which supplies the
  soft contributions with an effective cutoff.
  (3) The flavor symmetry breaking effects,
  and especially the distinction between different the twist DAs, are
  highlighted, compared with the nearly symmetric distribution
  Eq.(<ref>).

  
  
  
  
  
  

  
  

 §.§ Decay amplitudes

  
  The Feynman diagrams for the B^∗0 → D^+D^- decay
  with the pQCD approach are shown in Fig.<ref>, including the
  factorizable emission topologies (a,b) where one gluon links with
  the initial and the recoiled states, the
  nonfactorizable emission topologies (c,d) where one gluon is exchanged
  between the spectator quark and the emitted states, the
  factorizable annihilation topologies (e,f,i,j) where one gluon
  is conjoined with the final states, and the nonfactorizable annihilation
  topologies (g,h,k,l) where one gluon is transmitted between the
  initial and the final states.

  Generally, the amplitudes for the factorizable emission
  topologies in Fig.<ref>(a,b) can be written as the D
  meson decay constant and the space-like B^∗ → D
  transition form factor, and the amplitudes for the factorizable
  annihilation topologies in Fig.<ref>(e,f,i,j) can be
  written as the B^∗ meson decay constant and the time-like
  transition form factor between two charmed mesons.
  The amplitudes for the nonfactorizable topologies have more
  complicated structures, and can be written as the convolution
  integral of all participating meson wave functions.
  Compared with the contributions of the emission topologies
  in Fig.<ref>(a-d), the contributions of the annihilation
  topologies in Fig.<ref>(e-l) are assumed to be power
  suppressed, as stated in Ref.<cit.>.
  In addition, different topologies have different scales.
  The gluons of the emission topologies in Fig.<ref>(a-d)
  are time-like, while the gluons of the annihilation topologies
  in Fig.<ref>(e-l) are space-like.
  The gluon virtuality of creating a pair of heavy charm quarks
  from the vacuum for the annihilation topologies in Fig.<ref>(i-l),
  k_g^2 ≥ (2m_c)^2, should be much larger than
  that of producing a pair of light quarks for the annihilation
  topologies in Fig.<ref>(e-h).
  Thus, it is not hard to figure out that the contributions of the
  annihilation topologies in Fig.<ref>(i-l) might be very
  small relative to the others, because of the nature of the
  asymptotic freedom of the QCD at the unltrahigh energy.

  After a straightforward calculation, the amplitudes for the B_q^∗
  → DD decays are expressed as below.
  
   
    A(B^∗+→D^0D_q^+)
       
       =    F {
       V_cb^∗ V_cq [ a_1  A_a+b^LL+C_2  A_c+d^LL]
      +V_ub^∗ V_uq [ a_1  A_i+j^LL+C_2  A_k+l^LL]
       
            
      -V_tb^∗ V_tq [ (a_4+a_10)  A_a+b+i+j^LL
            +(a_6+a_8)  A_a+b^SP
                  
          + (C_3+C_9)  A_c+d+k+l^LL
            +(C_5+C_7)  A_c+d+k+l^SP] },

  
  
   
    A(B_q^∗0→D^0D^0)
       
       =    F {
       V_cb^∗ V_cq [ a_2  A_e+f^LL+C_1  A_g+h^LL]
      +V_ub^∗ V_uq [ a_2  A_i+j^LL+C_1  A_k+l^LL]
       
            
      -V_tb^∗ V_tq [ (a_3+a_9)  A_e+f+i+j^LL
            +(a_5+a_7)  A_e+f+i+j^LR
                  
          + (C_4+C_10)  A_g+h+k+l^LL
            +(C_6+C_8 )  A_g+h+k+l^LR] },

  
  
   
    A(B_q^∗0→D_q^-D_q^+)
       
       =    F {
       V_cb^∗ V_cq [ a_1  A_a+b^LL+C_2  A_c+d^LL
           + a_2  A_e+f^LL+C_1  A_g+h^LL]
       
            
      -V_tb^∗ V_tq [ (a_4+a_10)  A_a+b^LL
           + (a_6+a_8)  A_a+b^SP
           + (a_3+a_9)  A_e+f^LL
            
          + (C_3+C_9)  A_c+d^LL
          +  (C_5+C_7)  A_c+d^SP
          + (C_4+C_10)  A_g+h^LL
            
          + (a_5-1/2 a_7)  A_i+j^LR
          +  (C_5-1/2 C_7)  A_k+l^SP
          +  (C_6-1/2 C_8)  A_k+l^LR
            
          + (a_5+a_7)  A_e+f^LR
          +  (a_3+a_4-1/2 a_9-1/2 a_10)  A_i+j^LL
            
          + (C_6+C_8)  A_g+h^LR
          +  (C_3+C_4-1/2 C_9-1/2 C_10)  A_k+l^LL] },

  
  
   
    A(B_d^∗0→D^-D_s^+)
       
       =    F {
       V_cb^∗ V_cs [ a_1  A_a+b^LL+C_2  A_c+d^LL]
      -V_tb^∗ V_ts [ (a_4+a_10)  A_a+b^LL
            
          + (a_4-1/2 a_10)  A_i+j^LL
          +  (C_3-1/2 C_9)  A_k+l^LL
          +  (C_5-1/2 C_7)  A_k+l^SP
            
          + (a_6+a_8 )  A_a+b^SP
          +  (C_3+C_9 )  A_c+d^LL
          +  (C_5+C_7 )  A_c+d^SP] },

  
  
   
    A(B_s^∗0→D_s^-D^+)
       
       =    F {
       V_cb^∗ V_cd [ a_1  A_a+b^LL+C_2  A_c+d^LL]
      -V_tb^∗ V_td [ (a_4+a_10)  A_a+b^LL
            
          + (a_4-1/2 a_10)  A_i+j^LL
          +  (C_3-1/2 C_9)  A_k+l^LL
          +  (C_5-1/2 C_7)  A_k+l^SP
            
          + (a_6+a_8 )  A_a+b^SP
          +  (C_3+C_9 )  A_c+d^LL
          +  (C_5+C_7 )  A_c+d^SP] },

  
  
   
    A(B_d^∗0→D_s^-D_s^+)
       
       =    F {
       V_cb^∗ V_cd [ a_2  A_e+f^LL + C_1  A_g+h^LL]
       -V_tb^∗ V_td [ (a_3+a_9)  A_e+f^LL
            
          + (a_5+a_7 )  A_e+f^LR
          +  (C_4+C_10)  A_g+h^LL
          +  (C_6+C_8 )  A_g+h^LR
            
          + (a_3-1/2 a_9 )  A_i+j^LL
          +  (C_4-1/2 C_10)  A_k+l^LL
            
          + (a_5-1/2 a_7 )  A_i+j^LR
          +  (C_6-1/2 C_8 )  A_k+l^LR] },

  
  
   
    A(B_s^∗0→D^-D^+)
       
       =    F {
       V_cb^∗ V_cs [ a_2  A_e+f^LL + C_1  A_g+h^LL]
       -V_tb^∗ V_ts [ (a_3+a_9)  A_e+f^LL
            
          + (a_5+a_7 )  A_e+f^LR
          +  (C_4+C_10)  A_g+h^LL
          +  (C_6+C_8 )  A_g+h^LR
            
          + (a_3-1/2 a_9 )  A_i+j^LL
          +  (C_4-1/2 C_10)  A_k+l^LL
            
          + (a_5-1/2 a_7 )  A_i+j^LR
          +  (C_6-1/2 C_8 )  A_k+l^LR] },

  
   
    F  = √(2) G_F π C_F/N_c 
      f_B_q^∗  f_D̅  f_D  m_1 
      (ϵ_B_q^∗·p_D̅)
       ,

  
  
  where C_i is the Wilson coefficient; the parameter a_i is defined as
  
  
    a_i = {[  C_i+C_i+1/N_c  for odd i;; C_i+C_i-1/N_c  for even i, ].

  
  and A_m_1+m_2+⋯^n is
  an abbreviation for A_m_1^n +
  A_m_2^n + ⋯, where the
  subscript m_i corresponds to one of indices of
  Fig.<ref>; the superscript n refers
  to three possible Dirac structures Γ_1⊗Γ_2
  of the operators (q̅_1q_2)_Γ_1(q̅_3q_4)_Γ_2,
  namely n = LL for (V-A)⊗(V-A), n = LR
  for (V-A)⊗(V+A), and n = SP for -2(S-P)⊗(S+P).
  The explicit expressions of the building blocks A_m_i^n
  are collected in Appendix.

  
  

§ NUMERICAL RESULTS AND DISCUSSION

  
  In the rest frame of the B_q^∗ meson, the branching ratio is
  defined as
  
   
    Br = 1/24π p/m_B^∗^2Γ_B_q^∗ | A(B^∗→DD)|^2,

  
  where Γ_B_q^∗ is the full decay width of
   the B_q^∗ meson.

  
   
  

  The numerical values of some input parameters are listed in
  Table <ref>, where if it is not specified explicitly,
  their central values will be fixed as the default inputs.
  Besides, the full decay width of the B_q^∗ meson,
  Γ_B_q^∗, is also an essential parameter.
  Unfortunately, an experimental measurement on Γ_B_q^∗
  is unavailable now, because the soft photon from the
  B_q^∗ → B_qγ process is usually beyond
  the detection capability of electromagnetic calorimeters
  sitting at existing high energy colliders.
  It is well known that the electromagnetic radiation process
  B_q^∗ → B_qγ dominates the decay of
  the B_q^∗ meson.
  So, for the time being, the full decay width will be approximated
  by the radiative partial width, i.e., Γ_B_q^∗
  ≃ Γ(B_q^∗→B_qγ).
  At present, the information on Γ(B_q^∗→B_qγ)
  comes mainly from theoretical estimation.
  Theoretically, the partial decay width of the M1 transition (spin-flip)
  process has the expression <cit.>
  
   
    Γ(B_q^∗→B_qγ) =4/3 α  k_γ^3 μ^2_h,

  
  where α is the fine structure constant;
  k_γ = (m_B_q^∗^2-m_B_q^2)/2m_B_q^∗
  is the photon momentum in the rest frame of the B_q^∗ meson;
  μ_h is the M1 moment of the B_q^∗ meson.
  There are plenty of theoretical predictions on
  Γ(B_q^∗→B_qγ), for example,
  the numbers in Table 7 in Ref.<cit.> and
  Tables 3 and 4 in  Ref.<cit.>, but these
  estimation suffer from large uncertainties due to our
  insufficient understanding on the M1 moments of mesons.
  In principle, the M1 moment of a meson should be a
  combination of the M1 moments of the constituent quark and antiquark.
  For a heavy-light meson, the M1 moment of a heavy quark might be
  negligible relative to the M1 moment of a light quark, because
  it is widely assumed that the mass of a heavy quark is usually
  much larger than the mass of a light quark, and that the M1 moment
  is inversely proportional to the mass of a charged particle.
  With the M1 moment relations among light u, d, s quarks,
  |μ_u| > |μ_d| >
  |μ_s| <cit.>, one could expect to have
  Γ(B_u^∗→B_uγ) >
  Γ(B_d^∗→B_dγ) >
  Γ(B_s^∗→B_sγ),
  and so Γ_B_u^∗ > Γ_B_d^∗
  > Γ_B_s^∗.
  Of course, more details about the width Γ_B_q^∗
  is beyond the scope of this paper.
  In our calculation, in order to give a quantitative estimation of
  the branching ratios for the B^∗ → DD decays,
  we will fix
  
   
    Γ_B_u^∗ ∼Γ(B_u^∗→B_uγ) ∼   450 eV,

  
   
    Γ_B_d^∗ ∼Γ(B_d^∗→B_dγ) ∼   150 eV,

  
   
    Γ_B_s^∗ ∼Γ(B_s^∗→B_sγ) ∼   100 eV,

  
  which is basically consistent with the recent results
  in Refs.<cit.>.

  In order to investigate the effects from different DA models,
  we explore three scenarios,
  
  
  * Scenario I:
  ϕ_B^∗^v = Eq.(<ref>),
  ϕ_B^∗^t = Eq.(<ref>),
  ϕ_D̅,D^a = Eq.(<ref>)
  and ϕ_D̅,D^p = Eq.(<ref>).
  
  * Scenario II:
  ϕ_B^∗^v = ϕ_B^∗^t = Eq.(<ref>),
  and ϕ_D̅,D^a =
  ϕ_D̅,D^p = Eq.(<ref>).
  
  * Scenario III:
  ϕ_B^∗^v = ϕ_B^∗^t = Eq.(<ref>),
  ϕ_D̅,D^a =
  ϕ_D̅,D^p = Eq.(<ref>).
  

  Our numerical results are presented in Table <ref>,
  where the uncertainties come from the typical scale
  (1±0.1)t_i, mass m_c and m_b, and the CKM
  parameters, respectively.
  The following are some comments.

  
   
  

  (1)
  Generally, the B^∗ → DD decay modes may
  be divided into four classes. The tree contributions of classes
  A, B, C, D are proportional to the factors of
  V_cb^∗ V_cs a_1 ∼ Aλ^2 a_1,
  V_cb^∗ V_cd a_1 ∼ Aλ^3 a_1,
  V_cb^∗ V_cs a_2 ∼ Aλ^2 a_2,
  V_cb^∗ V_cd a_2 ∼ Aλ^3 a_2,
  respectively. Classes C and D are pure annihilation processes.
  There is a clear hierarchy of branching ratios, i.e.,
  Br(class A) ≳ 10^-9,
  Br(class B) ≳ 10^-11,
  Br(class C) ≳ 10^-13,
  Br(class D) ≲ 10^-13.

  In addition, for each class, the magnitude relations between the decay
  constants f_B_s^∗ > f_B_u,d^∗ and
  f_D_s > f_D_u,d, and relations among the decay widths
  Γ_B_s^∗ < Γ_B_d^∗
  < Γ_B_u^∗, result in the size relations
  among the branching ratios, i.e.,
  
   
    Br(B_s^∗0→D_s^-D_s^+)  >
       Br(B_d^∗0→D_d^-D_s^+)  >
       Br(B_u^∗+→D_u^0D_s^+)
       ,

  
   
    Br(B_s^∗0→D_s^-D_d^+)  >
       Br(B_d^∗0→D_d^-D_d^+)  >
       Br(B_u^∗+→D_u^0D_d^+)
       ,

  
   
    Br(B_d^∗0→D_s^-D_s^+)  >
       Br(B_d^∗0→D_u^0D_u^0)
       .

  

  (2)
  Due to the isospin symmetry, there are some approximate relations
  among the branching ratios, for example,
  
   
    Br(B_d^∗0→D_d^-D_s^+) / Br(B_u^∗+→D_u^0D_s^+)  ≈ Γ_B_u^∗/Γ_B_d^∗ ≈ 3
       ,

  
   
    Br(B_d^∗0→D_d^-D_d^+) / Br(B_u^∗+→D_u^0D_d^+)  ≈ Γ_B_u^∗/Γ_B_d^∗ ≈ 3
       ,

  
   
    Br(B_s^∗0→D_u^0D_u^0) / Br(B_s^∗0→D_d^-D_d^+)  ≈ 1
       .

  
  In addition, there are some other approximate relations,
  for example,
  
   
    Br(B_u^∗+→D_u^0D_s^+) / Br(B_u^∗+→D_u^0D_d^+)  ≈  f_D_s^2/ f_D_d^2 | V_cs|^2/| V_cd|^2,

  
   
    Br(B_d^∗0→D_d^-D_s^+) / Br(B_d^∗0→D_d^-D_d^+)  ≈  f_D_s^2/ f_D_d^2 | V_cs|^2/| V_cd|^2,

  
   
    Br(B_s^∗0→D_s^-D_s^+) / Br(B_s^∗0→D_s^-D_d^+)  ≈  f_D_s^2/ f_D_d^2 | V_cs|^2/| V_cd|^2.

  

  (3)
  Our study shows that
  (i) both the emission topologies and the annihilation topologies
  contribute to the decay channels of Classes A and B. Furthermore,
  the contributions
  from the emission topologies are dominant over those from the
  annihilation topologies.
  (ii) For the pure annihilation decay channels of Classes C and D,
  the factorizable contributions are color-suppressed, so the
  nonfactorizable contributions are the main ones.
  In addition, the interferences among different topologies are
  important.
  (iii) Compared with the contributions of the tree operators, the
  contributions of the penguin operators are small because of the
  suppression from the small Wilson coefficients.
  The contributions of the topologies in Fig.<ref>(i-l) are
  much less than those of the topologies in Fig.<ref>(e-h).
  (iv) For the decay channels with the final
  states D_qD_q, the contribution of the factorizable
  annihilation topology in Fig.<ref>(e) [Fig.<ref>(i)]
  is the same magnitude as that in Fig.<ref>(f) [Fig.<ref>(j)]
  due to the flavor symmetry.
  (v)
  The interferences between factorizable topologies in
  Fig.<ref>(e) and Fig.<ref>(f)
  [Fig.<ref>(i) and Fig.<ref>(j)] are constructive,
  while the interferences between nonfactorizable topologies in
  Fig.<ref>(g) and Fig.<ref>(h)
  [Fig.<ref>(k) and Fig.<ref>(l)] are destructive.

  
   
  

  (4)
  The B_s^∗0 → D_s^-D_s^+,
  B_d^∗0 → D_d^-D_s^+,
  B_u^∗+ → D_u^0D_s^+
  decays, belonging to class A, have relatively large branching ratios,
  Br(class A) ≳ 10^-9.
  The numbers of the B_q^∗ mesons in a data sample can be estimated by
  
   
    N(B_q^∗) =  L_ int × σ_bb̅ ×  B_q ×  B_B_q^∗/ B_q,

  
  
   
    B_B_q^∗ = 2 ×  B_B_q^∗B_q^∗
       +2 ×  B_B_q^∗B_q^∗π
       + B_B_q^∗B_q+ c.c.
       + B_B_q^∗B_qπ+ c.c
       +⋯,

  
  where L_ int is the integrated luminosity,
  σ_bb̅ denotes the bb̅ pair production cross section,
  B_q refers to the fragmentation fraction of (bq̅)(b̅q) events,
  and B_B_q^∗B_q^∗,
  B_B_q^∗B_q^∗X, ....
  represent the production fractions of specific modes (see Table <ref>).
  With a large production cross section of the process e^+e^- →
  bb̅ at the Υ(5S) peak σ_bb̅ =
  (0.340±0.016)  nb <cit.>
  and a high luminosity 8×10^35  cm^-2  s^-1
  at the forthcoming SuperKEKB <cit.>, it is expected that
  some 3.3×10^9 B_u,d^∗ and 1.2×10^9
  B_s^∗ mesons could be available per 10  ab^-1
  Υ(5S) dataset, corresponding to a few events of the
  B_d^∗0 → D_d^-D_s^+ and B_u^∗+
  → D_u^0D_s^+ decays and dozens of the
  B_s^∗0 → D_s^-D_s^+ decay.
  At high energy hadron colliders, for example, with a visible
  b-hadron cross section at the LHCb σ_bb̅
  ∼ 100 μb <cit.>,
  a similar ratio B_q at Tevatron and a similar ratio
  B_B_q^∗/ B_q at Υ(5S),
  some 9.2×10^13 B_u,d^∗
  and 1.9×10^13 B_s^∗ events per ab^-1
  dataset should be available at the LHCb, corresponding to more than
  10^5 (class A) B_s^∗0 → D_s^-D_s^+,
  B_d^∗0 → D_d^-D_s^+, and
  B_u^∗+ → D_u^0D_s^+
  decay events and over 10^4 (class B)
  B_s^∗0 → D_s^-D_d^+,
  B_d^∗0 → D_d^-D_d^+
  B_u^∗+ → D_u^0D_d^+ decay events,
  which are measurable at the future LHCb experiments.
  However, the search for pure annihilation processes (classes
  C and D) at LHC and SuperKEKB should still be very challenging.

  (5)
  Compared with the branching ratios ≳ 10^-5 for the B
  → DD decays <cit.>, the branching
  ratios for the B^∗ → DD decays are smaller by
  at least three orders of magnitude. This fact might imply that
  the background from the B^∗ → DD decays
  could be safely neglected when one analyzes the B →
  DD decays, but not vice versa, i.e., one of main
  pollution backgrounds for the B^∗ → DD
  decays would come from the B → DD decays,
  even if the invariant mass of the DD meson pair
  could be used to distinguish the B^∗ meson from the B meson
  experimentally.

  (6)
  For the B_q^∗ decays of classes A and B, our estimation of
  the branching ratios agrees well with that based on the naive
  factorization approach <cit.>.
  One of the important reasons is that these processes are all
  a_1-dominated (color favored), and in general, insensitive
  to nonfactorizable corrections to the HTME.
  Of course, one fact is clear that there are many theoretical
  uncertainties, especially, regarding the discrepancy among different DA
  scenarios, which results from our uncertain knowledge of
  the long-distance QCD effects and the underlying dynamics
  of low energy hadron interactions.
  Moreover, as aforementioned, there are large uncertainties
  of the decay width Γ_B_q^∗. With a different value
  of Γ_B_q^∗, the branching ratios
  in Table <ref> should be multiplied by the factors of
  450  eV/Γ_B_u^∗,
  150  eV/Γ_B_d^∗,
  100  eV/Γ_B_s^∗
  for the B_u^∗, B_d^∗, B_s^∗
  decays, respectively.
  In addition, many other factors, such as the final state
  interactions, models for the B^∗ and D meson wave functions[2],
  [2]In principle, one can do a global fit on the B^∗
  and D meson wave functions with experimental measurements in the future,
  analogous to that with the χ^2 method in Ref.<cit.>.
  The fitting will be a very time-consuming work,
  because the amplitudes for the B^∗ → DD decays
  are expressed as the multidimensional integral with the pQCD approach.
  In addition, there is no measurement report on the B^∗ →
  DD decays at the moment.
  higher order corrections to the HTME, and so on, are not
  carefully scrutinized here, but deserve much dedicated
  study. Our estimation may be just an order of magnitude.

  
  

§ SUMMARY

  
  With the running LHC and the forthcoming SuperKEKB, a large amount of
  B^∗ data should be in stock soon, which will make it seemingly
  possible to explore the B^∗ weak decays experimentally.
  A theoretical study is necessary in order to offer
  a timely reference, and is helpful in clearing up some of
  puzzles surrounding heavy meson weak decays.
  In this paper, we investigated the B^∗ →
  DD decays with the phenomenological pQCD approach.
  It is found that the B_s^∗0 → D_s^-D_s^+,
  B_d^∗0 → D_d^-D_s^+, and
  B_u^∗+ → D_u^0D_s^+
  decays have branching ratios ≳ 10^-9, and
  will be promisingly accessible at the future high luminosity
  experiments, with help of a sophisticated experimental analytical
  technique to effectively suppress or exclude the background
  from the B → DD decays.

  
  

§ ACKNOWLEDGMENTS

  The work is supported by the National Natural Science Foundation
  of China (Grant Nos. U1632109, 11547014 and 11475055).

  
  
  

§ AMPLITUDE BUILDING BLOCKS FOR B^∗ → DD DECAYS

  
  
  
   
    A_a^LL   =   ∫_0^1dx_1∫_0^1dx_2∫_0^∞b_1db_1∫_0^∞b_2db_2 α_s(t_a) 
      H_ef(α_e,β_a,b_1,b_2)
       
          ×  E_ef(t_a) ϕ_B^∗^v(x_1) {ϕ_D̅^a(x_2) 
       ( m_1^2 x̅_2+m_3^2 x_2 )
       + ϕ_D̅^p(x_2)  m_2 m_b},

  
   
    A_a^SP   =   
      -2 m_3 ∫_0^1dx_1∫_0^1dx_2∫_0^∞b_1db_1∫_0^∞b_2db_2 
      H_ef(α_e,β_a,b_1,b_2)
       
       ×   α_s(t_a)  E_ef(t_a) ϕ_B^∗^v(x_1) {ϕ_D̅^a(x_2) m_b
       + ϕ_D̅^p(x_2)  m_2 x_2},

  
  
  
   
    A_b^LL   =   ∫_0^1dx_1∫_0^1dx_2∫_0^∞b_1db_1∫_0^∞b_2db_2 α_s(t_b) 
      H_ef(α_e,β_a,b_2,b_1)
       
       ×   
      E_ef(t_b) {ϕ_B^∗^t(x_1) [
      ϕ_D̅^p(x_2)  2 m_1 m_2 x̅_1
     -ϕ_D̅^a(x_2) m_1 m_c]
       
       +   ϕ_B^∗^v(x_1) [
      ϕ_D̅^p(x_2)  2 m_2 m_c
     -ϕ_D̅^a(x_2) 
      ( m_2^2 x̅_1+m_3^2 x_1 ) ] },

  
  
   
    A_b^SP   =   
      2 m_3 ∫_0^1dx_1∫_0^1dx_2∫_0^∞b_1db_1∫_0^∞b_2db_2 α_s(t_b) 
      H_ef(α_e,β_a,b_2,b_1)
       
          × 
      E_ef(t_b) {ϕ_B^∗^v(x_1) [
      ϕ_D̅^a(x_2)  m_c
     -ϕ_D̅^p(x_2)  2 m_2]
     +ϕ_B^∗^t(x_1) ϕ_D̅^a(x_2) m_1 x_1},

  
  
  
   
    A_c^LL   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞db_1∫_0^∞b_2db_2∫_0^∞b_3db_3 α_s(t_c) 
      E_n(t_c)
       
       ×   
      H_en(α_e,β_c,b_3,b_2,b_1) {ϕ_B^∗^t(x_1) ϕ_D̅^p(x_2) ϕ_D^a(x_3)  m_1  m_2 (x_1-x_2)
       
           +
      ϕ_B^∗^v(x_1) ϕ_D̅^a(x_2) [
      ϕ_D^a(x_3) s (x_2-x̅_3)
     -ϕ_D^p(x_3)  m_3  m_c] },

  
  
   
    A_c^SP   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞db_1∫_0^∞b_2db_2∫_0^∞b_3db_3 
      H_en(α_e,β_c,b_3,b_2,b_1)
       
       ×   
       E_n(t_c) α_s(t_c) {ϕ_B^∗^v(x_1) ϕ_D̅^p(x_2)  m_2 [
      ϕ_D^p(x_3)   m_3 (x_2-x̅_3)
     -ϕ_D^a(x_3)   m_c]
       
           +
      ϕ_B^∗^t(x_1) ϕ_D̅^a(x_2)  m_1 [
      ϕ_D^a(x_3)  m_c
     +ϕ_D^p(x_3)  m_3  (x̅_3-x_1)
       ] },

  
  
  
   
    A_d^LL   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞db_1∫_0^∞b_2db_2∫_0^∞b_3db_3 α_s(t_d) 
      E_n(t_d)
       
       ×   
      H_en(α_e,β_d,b_3,b_2,b_1) ϕ_D^a(x_3) {ϕ_B^∗^t(x_1) ϕ_D̅^p(x_2)  m_1 m_2 (x_1-x_2)
       
           +
       ϕ_B^∗^v(x_1) ϕ_D̅^a(x_2) 
       (2 m_2^2 x_2+s x_3-t x_1) },

  
  
   
    A_d^SP   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞db_1∫_0^∞b_2db_2∫_0^∞b_3db_3 α_s(t_d) 
      E_n(t_d)
       
       ×   
      H_en(α_e,β_d,b_3,b_2,b_1) ϕ_D^p(x_3) {ϕ_B^∗^v(x_1) ϕ_D̅^p(x_2)  m_2 m_3 (x_3-x_2)
       
           +
       ϕ_B^∗^t(x_1) ϕ_D̅^a(x_2) 
       m_1 m_3 (x_1-x_3) },

  
  
  
   
    A_e^LL   =    A_e^LR =∫_0^1dx_2∫_0^1dx_3∫_0^∞b_2db_2∫_0^∞b_3db_3 
      H_af(α_q,β_e,b_2,b_3)
       
       ×   α_s(t_e)  E_af(t_e) {ϕ_D^p(x_3) 2 m_3 [
      ϕ_D̅^a(x_2) m_c
     +ϕ_D̅^p(x_2) m_2 x̅_2]
       
          
     -ϕ_D^a(x_3) [
      ϕ_D̅^a(x_2) (m_1^2 x_2+m_3^2x̅_2 )
     +ϕ_D̅^p(x_2) m_2 m_c] },

  
  
  
   
    A_f^LL   =    A_f^LR =∫_0^1dx_2∫_0^1dx_3∫_0^∞b_2db_2∫_0^∞b_3db_3 
      H_af(α_q,β_f,b_3,b_2)
       
       ×   α_s(t_f)  E_af(t_f) {ϕ_D̅^p(x_2) 2 m_2 [
      ϕ_D^a(x_3) m_c
     +ϕ_D^p(x_3) m_3 x̅_3]
       
          
     -ϕ_D̅^a(x_2) [
      ϕ_D^a(x_3) (m_1^2 x_3+m_2^2x̅_3 )
     +ϕ_D^p(x_3) m_3 m_c] },

  
  
  
   
    A_g^LL   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3 
      H_an(α_q,β_g,b_1,b_2,b_3)
       
       ×   
       E_n(t_g) α_s(t_g) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) [ ϕ_B^∗^v(x_1) 
      (s x_2+2 m_3^2 x_3-u x̅_1)
       
          
      +ϕ_B^∗^t(x_1) m_1 m_b]
      +ϕ_B^∗^v(x_1) ϕ_D̅^p(x_2) ϕ_D^p(x_3)  m_2 m_3 (x_2-x_3) },

  
  
   
    A_g^LR   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3 
      H_an(α_q,β_g,b_1,b_2,b_3)
       
       ×   
       E_n(t_g) α_s(t_g) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) [ ϕ_B^∗^v(x_1) 
      (t x̅_1-2 m_2^2 x_2-s x_3)
       
          
      -ϕ_B^∗^t(x_1) m_1 m_b]
      +ϕ_B^∗^v(x_1) ϕ_D̅^p(x_2) ϕ_D^p(x_3)  m_2 m_3 (x_2-x_3) },

  
  
  
   
    A_h^LL   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3 
      E_n(t_h)
       
       ×   α_s(t_h) ϕ_B^∗^v(x_1) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) 
      (2 m_2^2 x_2+s x_3-t x_1)
       
       +   ϕ_D̅^p(x_2) ϕ_D^p(x_3) 
      m_2 m_3 (x_3-x_2) } 
      H_an(α_q,β_h,b_1,b_2,b_3)
       ,

  
  
   
    A_h^LR   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3 
      E_n(t_h)
       
       ×   α_s(t_h) ϕ_B^∗^v(x_1) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) 
      (u x_1-s x_2- 2 m_3^2 x_3)
       
       +   ϕ_D̅^p(x_2) ϕ_D^p(x_3) 
      m_2 m_3 (x_3-x_2) } 
      H_an(α_q,β_h,b_1,b_2,b_3)
       ,

  
  
  
   
    A_i^LL   =    A_i^LR =∫_0^1dx_2∫_0^1dx_3∫_0^∞b_2db_2∫_0^∞b_3db_3 α_s(t_i) 
      H_af(α_c,β_i,b_2,b_3)
       
          ×  E_af(t_i) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) 
      (m_1^2 x̅_2+m_3^2 x_2)
     -ϕ_D̅^p(x_2) ϕ_D^p(x_3) 
      2 m_2 m_3 x_2},

  
  
  
   
    A_j^LL   =    A_j^LR =∫_0^1dx_2∫_0^1dx_3∫_0^∞b_2db_2∫_0^∞b_3db_3 α_s(t_j) 
      H_af(α_c,β_j,b_3,b_2)
       
          ×  E_af(t_j) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) 
      (m_1^2 x̅_3+m_2^2 x_3)
     -ϕ_D̅^p(x_2) ϕ_D^p(x_3) 
      2  m_2 m_3 x_3},

  
  
  
   
    A_k^LL   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3 
      H_an(α_c,β_k,b_1,b_2,b_3)
       
       ×   
       E_n(t_k) α_s(t_k) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) [ ϕ_B^∗^v(x_1) 
      (2 m_2^2 x_2+s x_3-t x_1)
       
          
      -ϕ_B^∗^t(x_1) m_1 m_b]
      +ϕ_B^∗^v(x_1) ϕ_D̅^p(x_2) ϕ_D^p(x_3)  m_2 m_3 (x_3-x_2) },

  
  
   
    A_k^LR   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3 
      H_an(α_c,β_k,b_1,b_2,b_3)
       
       ×   
       E_n(t_k) α_s(t_k) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) [ ϕ_B^∗^v(x_1) 
      (u x_1-s x_2-2 m_3^2 x_3)
       
          
      +ϕ_B^∗^t(x_1) m_1 m_b]
      +ϕ_B^∗^v(x_1) ϕ_D̅^p(x_2) ϕ_D^p(x_3)  m_2 m_3 (x_3-x_2) },

  
  
   
    A_k^SP   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3 
      H_an(α_c,β_k,b_1,b_2,b_3)
       
       ×   
       E_n(t_k) {ϕ_D̅^a(x_2) ϕ_D^p(x_3) 
       m_3 [ ϕ_B^∗^v(x_1)  m_b
       +ϕ_B^∗^t(x_1) m_1 (x_1-x_3) ]
       
          
      +ϕ_D̅^p(x_2) ϕ_D^a(x_3) 
       m_2 [ ϕ_B^∗^v(x_1) m_b
      +ϕ_B^∗^t(x_1) m_1 (x_1-x_2)
        ] } α_s(t_k)
       ,

  
  
  
   
    A_l^LL   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3 
      E_n(t_l)
       
       ×   α_s(t_l) ϕ_B^∗^v(x_1) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) 
      (u x_1-s x̅_2-2 m_3^2 x̅_3)
       
       +   ϕ_D̅^p(x_2) ϕ_D^p(x_3) 
      m_2 m_3 (x_2-x_3) } 
      H_an(α_c,β_l,b_1,b_2,b_3)
       ,

  
  
   
    A_l^LR   =   1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3 
      E_n(t_l)
       
       ×   α_s(t_l) ϕ_B^∗^v(x_1) {ϕ_D̅^a(x_2) ϕ_D^a(x_3) 
      (2 m_2^2 x̅_2+s x̅_3-t x_1)
       
       +   ϕ_D̅^p(x_2) ϕ_D^p(x_3) 
      m_2 m_3 (x_2-x_3) } 
      H_an(α_c,β_l,b_1,b_2,b_3)
       ,

  
  
   
    A_l^SP   =   m_1/N_c ∫_0^1dx_1∫_0^1dx_2∫_0^1dx_3∫_0^∞b_1db_1∫_0^∞b_2db_2∫_0^∞db_3
       ×    E_n(t_l) α_s(t_l) ϕ_B^∗^t(x_1) {ϕ_D̅^a(x_2) ϕ_D^p(x_3) 
       m_3  (x̅_3-x_1)
       
       +   ϕ_D̅^p(x_2) ϕ_D^a(x_3) 
       m_2 (x̅_2-x_1) } 
      H_an(α_c,β_l,b_1,b_2,b_3)
       ,

  
  
  where the subscript i of A_i^j corresponds to
  the indices of Fig.<ref>; the superscript j refers
  to three possible Dirac structures Γ_1⊗Γ_2
  of the operators (q̅_1q_2)_Γ_1(q̅_3q_4)_Γ_2,
  namely j = LL for (V-A)⊗(V-A), j = LR
  for (V-A)⊗(V+A), and j = SP for -2(S-P)⊗(S+P).

  The function H_i and the Sudakov factor E_i are defined as
  
  
   
    H_ef(α_e,β,b_m,b_n)  = 
       K_0(b_m√(-α_e)) {θ(b_m-b_n) K_0(b_m√(-β)) 
       I_0(b_n√(-β))
       + (b_m↔ b_n) },

  
  
   
    H_en(α_e,β,b_3,b_2,b_1)    =   {θ(-β)  K_0(b_3√(-β))
      +π/2 θ(β) [ i J_0(b_3√(β))
       - Y_0(b_3√(β)) ] }
          × {θ(b_2-b_3)  K_0(b_2√(-α_e)) 
       I_0(b_3√(-α_e)) + (b_2↔ b_3)
       } δ(b_1-b_2)
       ,

  
  
   
    H_af(α,β,b_m,b_n)
          =   {θ(b_m-b_n) [ i J_0(b_m√(β))
       - Y_0(b_m√(β)) ]
       J_0(b_n√(β))
       
          
      +  (b_m↔ b_n) } { i J_0(b_n√(α))
       - Y_0(b_n√(α)) } π^2/4,

  
  
   
    H_an(α,β,b_1,b_2,b_3)
          =   π/2 {θ(-β)  K_0(b_1√(-β))
      + π/2θ(β) [ i J_0(b_1√(β))
       - Y_0(b_1√(β)) ] }
          × {θ(b_1-b_2) [ i J_0(b_1√(α))
       - Y_0(b_1√(α)) ]
       J_0(b_2√(α))
       + (b_1↔ b_2) } δ(b_2-b_3)
       ,

  
   
    E_ef(t) = exp{ -S_B^∗(t)-S_D̅(t) },

  
  
   
    E_af(t) = exp{ -S_D̅(t)-S_D(t) },

  
  
   
    E_n(t) = exp{ -S_B^∗(t)-S_D̅(t)-S_D(t) },

  
   
    S_B^∗(t)  =  s(x_1,b_1,p_1^+)
      +2∫_1/b_1^tdμ/μγ_q,

  
  
    S_D̅(t)  =  s(x_2,b_2,p_2^+) + s(x̅_2,b_2,p_2^+)
      +2∫_1/b_2^tdμ/μγ_q,

  
  
    S_D(t)  =  s(x_3,b_3,p_3^+) + s(x̅_3,b_3,p_3^+)
      +2∫_1/b_3^tdμ/μγ_q,

  
  
  where the subscript i = ef, en, af, an corresponds
  to the factorizable emission topologies, the nonfactorizable emission
  topologies, the factorizable annihilation topologies, and the nonfactorizable
  annihilation topologies, respectively;
  I_0, J_0, K_0 and Y_0 are Bessel functions;
  γ_q = -α_s/π is the quark anomalous
  dimension; the expression of s(x,b,Q) can be found in of
  Ref.<cit.>;
  α and β are the virtualities of gluon and quarks.
  the subscript of β_i corresponds to the indices of
  Fig.<ref>.
  The definitions of the particle virtuality and typical
  scale t_i are given as follows.
  
  
   
    α_e = x_1^2 m_1^2+x_2^2 m_2^2-x_1 x_2 t
       ,

  
  
   
    α_q = x_2^2  m_2^2+x_3^2 m_3^2+x_2 x_3 s
       ,

  
  
   
    α_c = x̅_2^2 m_2^2+x̅_3^2 m_3^2
      +x̅_2 x̅_3 s
       ,

  
  
  
   
    β_a = m_1^2 + x_2^2 m_2^2-x_2 t-m_b^2,

  
  
   
    β_b = m_2^2+x_1^2 m_1^2-x_1 t - m_c^2,

  
  
   
    β_c = α_e + x̅_3^2 m_3^2
      -x_1 x̅_3 u +x_2 x̅_3 s -m_c^2,

  
  
   
    β_d = α_e+x_3^2 m_3^2-x_1 x_3 u+x_2 x_3 s
       ,

  
  
  
   
    β_e = m_3^2+x_2^2 m_2^2+x_2 s-m_c^2 > 0
       ,

  
  
   
    β_f = m_2^2+x_3^2 m_3^2+x_3 s-m_c^2 > 0
       ,

  
  
   
    β_g = α_q+x̅_1^2 m_1^2
                   - x̅_1 x_2 t - x̅_1 x_3 u -m_b^2,

  
  
   
    β_h = α_q+x_1^2 m_1^2
                   - x_1 x_2 t - x_1 x_3 u
       ,

  
  
  
   
    β_i = m_3^2+x̅_2^2 m_2^2+x̅_2 s > 0
       ,

  
  
   
    β_j = m_2^2+x̅_3^2 m_3^2+x̅_3 s > 0
       ,

  
  
   
    β_k = α_c+x̅_1^2 m_1^2
                   - x̅_1 x̅_2 t
                   - x̅_1 x̅_3 u -m_b^2,

  
  
   
    β_l = α_c+x_1^2 m_1^2
                   - x_1 x̅_2 t - x_1 x̅_3 u
       ,

  
  
  
  
   
    t_a,b = max{√(-α_e),√(|β_a,b|),1/b_1,1/b_2},

  
  
   
    t_c,d = max{√(-α_e),√(|β_c,d|),1/b_2,1/b_3},

  
  
   
    t_e,f = max{√(α_q),√(|β_e,f|),1/b_2,1/b_3},

  
  
   
    t_g,h = max{√(α_q),√(|β_g,h|),1/b_1,1/b_2},

  
  
   
    t_i,j = max{√(α_c),√(|β_i,j|),1/b_2,1/b_3},

  
  
   
    t_k,l = max{√(α_c),√(|β_k,l|),1/b_1,1/b_2}.

  
  
  

  
  99
  pdg
          C. Patrignani et al. (Particle Data Group), Chin. Phys. C 40, 100001 (2016).
  1509.05049
          B. Grinstein, J. Camalich, Phys. Rev. Lett. 116, 141801 (2016).
  1601.03386
          G. Xu et al., Eur. Phys. J. C 76, 583 (2016).
  1605.01630
          Q. Chang et al., Int. J. Mod. Phys. A 30, 1550162 (2015).
  1605.01629
          Q. Chang et al., Adv. in High Energy Phys. 2015, 767523 (2015).
  1605.01631
          Q. Chang et al., Eur. Phys. J. C 76, 523 (2016).
  1606.09071
          Q. Chang et al., Nucl. Phys. B 909, 921 (2016).
  zpc29
          M. Wirbel, B. Stech, M. Bauer, Z. Phys. C 29, 637 (1985).
  qcdf1
          M. Beneke et al., Phys. Rev. Lett. 83, 1914 (1999).
  qcdf2
          M. Beneke et al., Nucl. Phys. B 591, 313 (2000).
  qcdf3
          M. Beneke et al., Nucl. Phys. B 606, 245 (2001).
  plb.509.263
          D. Du, D. Yang, G. Zhu, Phys. Lett. B 509, 263 (2001).
  prd64.014036
           D. Du, D. Yang, G. Zhu, Phys. Rev. D 64, 014036 (2001).
  pqcd1
          H. Li, Phys. Rev. D 52, 3958 (1995).
  pqcd2
          C. Chang,  H. Li, Phys. Rev. D 55, 5577 (1997).
  pqcd3
          T. Yeh, H. Li, Phys. Rev. D 56, 1615 (1997).
  9512380
          G. Buchalla, A. Buras, M. Lautenbacher, Rev. Mod. Phys. 68, 1125, (1996).
  npb133
          D. Fakirov, B. Stech, Nucl. Phys. B 133, 315 (1978).
  zpc34
          M. Bauer, B. Stech, M. Wirbel, Z. Phys. C 34, 103 (1987).
  prd22
          G. Lepage, S. Brodsky, Phys. Rev. D 22, 2157 (1980).
  scet1
          C. Bauer et al.,  Phys. Rev. D 63, 114020 (2001).
  scet2
          C. Bauer, D. Pirjol, I. Stewart, Phys. Rev. D 65, 054022 (2002).
  scet3
          C. Bauer et al., Phys. Rev. D 66, 014017 (2002).
  scet4
          M. Beneke et al., Nucl. Phys. B 643, 431 (2002).
  npb529
          P. Ball, V. Braun, Y. Koike, K. Tanaka, Nucl. Phys. B 529, 323 (1998).
  prd65
          T. Kurimoto, H. Li, A. Sanda, Phys. Rev. D 65, 014007 (2001).
  plb752sun
          J. Sun et al., Phys. Lett. B 752, 322 (2016).
  ijmpa31yang
          Y. Yang et al., Int. J. Mod. Phys. A 31, 1650146 (2016).
  prd78.014018
          R. Li, C. Lü, H. Zou, Phys. Rev. D 78, 014018 (2008).
  prd81.034006
          R. Li et al., Phys. Rev. D 81, 034006 (2010).
  jpg31.273
          H. Zou et al., J. Phys. G 37, 015002 (2010).
  jpg37
          Y. Li, C. Lü, Z. Xiao, J. Phys. G 31, 273 (2005).
  prd91.114509
          B. Colquhoun et al. (HPQCD Collaboration), Phys. Rev. D 91, 114509 (2015).
  uds
          A. Kamal, Particle Physics, Springer, 2014, p. 297, p. 298.
  jhep1404.177
          C. Cheung, C. Hwang, JHEP, 1404, 177 (2014).
  epja52.90
          V. Šimonis, Eur Phys J A 52, 90 (2016).
  epjc74.3026
          Ed. A. Bevan et al., Eur. Phys. J. C 74, 3026 (2014). 
  1412.7515
          Y. Amhis et al. (Heavy Flavor Averaging Group), arXiv:1412.7515[hep-ex]. 
  1002.5012
          A. Akeroyd et al., arXiv:1002.5012[hep-ex].
  crp16.435
          T. Gershon, M. Needham, Comptes Rendus Physique 16, 435 (2015).
  
  