The Next-to-Leading Order Corrections to Top Quark Decays to Heavy Quarkonia
aa r X i v : . [ h e p - ph ] M a y The Next-to-Leading Order Corrections to Top Quark Decays toHeavy Quarkonia
Peng Sun , Li-Ping Sun , Cong-Feng Qiao , ∗ College of Physical Sciences, Graduate University of Chinese Academy of SciencesYuQuan Road 19A, Beijing 100049, China and Theoretical Physics Center for Science Facilities (TPCSF), CASYuQuan Road 19B, Beijing 100049, China
Abstract
The decay widths of top quark to S-wave b ¯ c and b ¯ b bound states are evaluated at the next-to-leading(NLO) accuracy in strong interaction. Numerical calculation shows that the NLO cor-rections to these processes are remarkable. The quantum chromodynamics(QCD) renormalizationscale dependence of the results is obviously depressed, and hence the uncertainties lying in theleading order calculation are reduced. PACS number(s): ∗ corresponding author . INTRODUCTION Since predicted by the Standard Model (SM)[1–3], top quark has become an importantrole in high energy physics due to its large mass, which is close to the electroweak symmetrybreaking scale [4]. A great deal of researches focusing on top quark physics have beenperformed after its discovery in 1995 in the Fermilab [5]. On the experiment aspect, withthe running of the Tevatron and forthcoming LHC, the lack of adequate events will not be anobstacle for the top quark physics study. According to Ref. [6], at the LHC 10 ∼ t ¯ t pairscan be obtained per year, so this enables people to measure various top quark decay channels.Meanwhile, the copious production of the top quarks supplies also a great number of bottomquark mesons since the dominant top quark decay channel is t → b + W + . Therefore,the bottom quark meson production in top quark decays may stand as an important andindependent means for the study of heavy meson physics and the test of perturbative QCD(pQCD).As the known heaviest mesons, bottomonia and B c ( b ¯ c or ¯ bc ) possess particular meaningin the study of heavy flavor physics. The LHCb as a detector specifically for the heavyflavor study at the LHC will supply copious B c and Υ data for this aim. Theoretically, thedirect hadroproduction of B c and Υ was studied in the literature [7–9]. In addition to the“direct” production, “indirect” process as in top quark decays may stand as an independentand important source for B c and Υ production. Since the top quark’s lifetime is too shortto form a bound state [10], the B c and Υ production involved scheme in top quark decaysis less affected by the non-pertubative effects than in other processes. In Ref. [11], the topquark decays into Υ and ¯ B ∗ c at the Born level was evaluated. Recently, the S- and P-wave B c meson productions in top quark decays were fully evaluated, including the color-octetcontributions, at the leading order accuracy of QCD by Chang et al. [12].Considering the importance of investigating B c and Υ in the study of perturbative Quan-tum Chromodynamics (pQCD) and potential model, it is reasonable and interesting to eval-uate the production rates of these mesons in top quark decays at the next-to-leading order(NLO) accuracy of pQCD. At the bottom quark and charm quark mass scales the strongcoupling is not very small, therefore the higher order corrections are usually large. On theother hand, in the processes of top quark decays into ¯ B c (Υ), the t → bc (¯ b ) + c ( b ) + W + ,2here exist large scale uncertainties in the tree level calculation [13]. The NLO correctionsshould in principle minimize it and give a more precise prediction. To calculate the ¯ B c andΥ production rates in top quark decays at the NLO accuracy are the aims of this work. Inour calculation, both of the S-wave spin-singlet and -triplet states are taken into account,i.e., ¯ B ∗ c , ¯ B c , Υ and η b . To deal with the non-perturbative effects, the non-relativistic QCD(NRQCD) [14] effective theory is employed. The calculation will be performed at the NLOin pQCD expansion, but at leading order in relativistic expansion, that is in the expansionof v , the relativistic velocity of heavy quarks inside bound states.The paper is organized as follows: after the Introduction, in section II we explain thecalculation of leading order decay width. In section III, virtual and real QCD correctionsto Born level result are evaluated. In section IV, the numerical calculation for concernedprocesses at NLO accuracy of pQCD is performed, and the scale dependence of the resultsis shown. The last section is remained for a brief summary and conclusions. II. CALCULATION OF THE BORN LEVEL DECAY WIDTH
At the leading order in α s , there are two Feynman Diagrams for each meson production,which are shown in Figure 1. For the convenience of analytical calculation, taking ¯ B c as anexample, the momentum of each particle is assigned as: p = p t , p = p b , p = p c , p = p c , p = p W + , p = p + p , p = m b m c p . For bottomonium, the only difference is that p and p represent the momenta of anti-bottom quark and bottom quark which are produced ingluon splitting.Of the ¯ B c and ¯ B ∗ c production in top quark decays, i.e. t ( p ) → ¯ B c / ¯ B ∗ c ( p ) + c ( p ) + W + ( p ) , (1)we employ the following commonly used projection operators for quarks hadronization: v ( p ) u ( p ) −→ √ iγ ( p + m b + m c ) × q m b + m c ψ ¯ B c (0) ⊗ (cid:18) c √ N c (cid:19) (2)and v ( p ) u ( p ) −→ √ ǫ ¯ B ∗ c ( p + m b + m c ) × q m b + m c ψ ¯ B ∗ c (0) ⊗ (cid:18) c √ N c (cid:19) . (3)3 B c ( ¯ B ∗ c ) cW + t t ¯ B c ( ¯ B ∗ c ) cW + ( a ) ( b ) FIG. 1: The leading order Feynman diagrams for ¯ B c and ¯ B ∗ c production in top quark decays. Here, ε ¯ B ∗ c is the polarization vector of ¯ B ∗ c with p · ε = 0, c stands for the unit colormatrix, and N c = 3 for QCD. The nonperturbative parameters ψ ¯ B c (0) and ψ ¯ B ∗ c (0) are theSchr¨odinger wave functions at the origin of b ¯ c bound states, and in the non-relativistic limit ψ ¯ B c (0) = ψ ¯ B ∗ c (0). In our calculation, the non-relativistic relation m ¯ B c = m ¯ B ∗ c = m b + m c isalso adopted.The LO amplitudes for ¯ B c production can then be readily obtained with above prepara-tions. They are: M a = πα s gψ ¯ B c (0) V tb C F δ j,k p m ¯ B c ¯ u ( p ) γ µ iγ ( p + m ¯ B c ) γ µ ( p + p + m b )( p + p ) − m b ǫ ( p )(1 − γ )( p + p ) u ( p ) , (4)and M b = πα s gψ ¯ B c (0) V tb C F δ j,k p m ¯ B c ¯ u ( p ) γ µ iγ ( p + m ¯ B c ) ǫ ( p )(1 − γ )( p + p ) ( p + p + m t )( p + p ) − m t γ µ u ( p ) . (5)Here, j , k are color indices, C F = 4 / SU (3) color structure. For ¯ B ∗ c production, the amplitudes can be obtained by simply substituting iγ ( p + m ¯ B c ) with ǫ ¯ B ∗ c ( p + m b + m c ) in above expressions.The Born amplitude of the processes shown in Fig.1 is then M Born = M a + M b , andsubsequently, the decay width at leading order reads:dΓ Born = 12 m t
12 1 N c X |M Born | dPS . (6)Here, P represents the sum over polarizations and colors of the initial and final particles, and N c come from spin and color average of initial t quark, dPS stands for the integrants4f three-body phase space, whose concrete form isdPS = 132 π m t d s d s , (7)where s = ( p + p ) = ( p − p ) and s = ( p + p ) = ( p − p ) are Mandelstam variables.The upper and lower bounds of the above integration are s max = q f [ m t , s , m B c ] · f [ s , m c , m W ] + [ m t − s − ( m b + m c ) ]( s + m c − m W )2 s + m B c + m c , (8) s min = − q f [ m t , s , m B c ] · f [ s , m c , m W ] − [ m t − s − ( m b + m c ) ]( s + m c − m W )2 s + m B c + m c (9)and s max = [ m t − ( m b + m c )] , s min = ( m c + m W ) (10)with f [ x, y, z ] = ( x − y − z ) − yz . (11) III. THE NEXT-TO-LEADING ORDER CORRECTIONS
At the next-to-leading order, the top quark decays to ¯ B c and Υ include the virtualand real QCD corrections to the leading order process, as shown in Figs.2-5. With virtualcorrections, the decay widths at the NLO can be formulated asdΓ V irtual = 12 m t
12 1 N c X M ∗ Born M V irtual )dPS . (12)The ultraviolet(UV) and infrared(IR) divergences usually exist in virtual corrections. Weuse the dimensional regularization scheme to regularize the UV and IR divergences, similaras performed in Ref.[15], and the Coulomb divergence is regularized by the relative velocity v . In dimensional regularization, γ is difficult to deal with. In this calculation, we adopt5 IG. 2: The self-energy diagrams in virtual corrections. the Naive scheme, that is, γ anticommutates with each γ µ matrix in d-dimension space-time, { γ , γ µ } = 0. The UV divergences exist merely in self-energy and triangle diagrams,which can be renormalized by counter terms. The renormalization constants include Z , Z , Z m , and Z g , corresponding to quark field, gluon field, quark mass, and strong couplingconstant α s , respectively. Here, in our calculation the Z g is defined in the modified-minimal-subtraction (MS) scheme, while for the other three the on-shell (OS) scheme is adopted,which tells δZ OSm = − C F α s π (cid:20) ǫ UV − γ E + ln 4 πµ m + 43 + O ( ǫ ) (cid:21) ,δZ OS = − C F α s π (cid:20) ǫ UV + 2 ǫ IR − γ E + 3 ln 4 πµ m + 4 + O ( ǫ ) (cid:21) ,δZ OS = α s π (cid:20) ( β − C A )( 1 ǫ UV − ǫ IR ) + O ( ǫ ) (cid:21) ,δZ MSg = − β α s π (cid:20) ǫ UV − γ E + ln 4 π + O ( ǫ ) (cid:21) . (13)Here, β = (11 / C A − (4 / T f n f is the one-loop coefficient of the QCD beta function; n f = 5 is the number of active quarks in our calculation; C A = 3 and T F = 1 / µ is the renormalization scale.In virtual corrections, IR divergences remain in the triangle and box diagrams. Ofall the triangle diagrams, only two have IR divergences, which are denoted as TriangleN7and TriangleN9 in Fig.3. Of the diagrams in Fig.4, BoxN3 has no IR singularity, whileBoxN4 and PentagonN9 have Coulomb singularities and PentagonN9 possesses ordinary6 IG. 3: The triangle diagrams in virtual corrections.
IR singularity as well. The remaining diagrams all have IR singularities, while the combi-nations BoxN2 + BoxN6, BoxN1 + PentagonN8 + TriangleN9, BoxN5 + TriangleN7 are IRfinite. The Coulomb singularities belonging to BoxN4 and PentagonN9 can be regularizedby the relative velocity v . After regularization procedure, the ǫ term will be canceled outby the counter terms of external quarks which form the ¯ B c or Υ, while the v term will bemapped onto the wave functions of the concerned heavy mesons. The remaining IR singu-larities in BoxN7 and BoxN9 are canceled by the corresponding parts in real corrections. Inthe end, the IR and Coulomb divergences in virtual corrections can be expressed asdΓ IR,Coulombvirtual = dΓ
Born α s π (cid:20) π v − ǫ − p t · p c x s ln x s m c m t (1 − x s ) 1 ǫ (cid:21) , (14)with p t = p , p c = p and x s = − √ − m c m t / ( m c m t − p t · p c )1+ √ − m c m t / ( m c m t − p t · p c ) . Here, in this work ǫ in factrepresents ǫ − γ E + ln(4 πµ ).Of our concerned processes, there are 12 different diagrams in real correction, as shown7 IG. 4: The box and pentagon diagrams in virtual corrections. in Fig.5. Among them, RealN2, RealN3, RealN8, and RealN9 are IR-finite, meanwhile thecombinations of RealN1 + RealN5 and RealN10 + RealN11 exibit no IR singularities as well,due to the reasons of gluon connecting to the b or ¯ c quark of final ¯ B c or Υ. The remainingdiagrams, RealN4, RealN6, RealN7, and RealN12 are not IR singularity free. To regularizethe IR divergence, we enforce a cut on the gluon momentum, the p . The gluon with energy p < δ is considered to be soft, while p > δ is thought to be hard. The δ is a small quantitywith energy-momentum unit. In this case, the IR term of the decay width can then bewritten as: dΓ IRReal = 12 m t
12 1 N c X |M Real | dPS | soft , (15)where dPS is the four-body phase space integrants for real correction. Under the conditionof p < δ , in the Eikonal approximation we obtaindPS | soft = dPS d p (2 π ) p | p <δ (16)8 IG. 5: The real correction Feynman diagrams that contribute to the production of ¯ B c or ¯ B ∗ c . In the small δ limit, the IR divergent terms in real correction can therefore be expressed asdΓ IRReal = dΓ
Born α s π (cid:26)(cid:18) ǫ − Log( δ ) (cid:19) (cid:20) p t · p c x s ln x s m c m t (1 − x s ) (cid:21) + finite terms (cid:27) . (17)Here, the Log( δ ) involved terms will be canceled out by the δ -dependent terms in the hardsector of real corrections. Referring to the Eq.(14), it is obvious that the IR divergentterms in real and virtual corrections cancel with each other. In case of hard gluons in realcorrection, the decay width readsdΓ hardReal = 12 m t
12 1 N c X |M Real | dPS | hard . (18)In this case the phase space dPS | hard can be written asdPS | hard = 2(4 π ) p ( sy + m c − m W ) − sym c y Z p p − d p Z − d cos θ c Z π d φ c × (Z p − δ d p Z y + y − d y + Z p p − d p Z y +( mc + mW )2 s d y ) (19)with p − = m b + m c , (20)9 = s + m b − m W + 2 m b m c − m W · m c √ s , (21) p − = s + m b − m W + 2 m b m c − m W · m c − √ sp √ s − p + 2 p |−→ p | , (22) p = s + m b − m W + 2 m b m c − m W · m c − √ sp √ s − p − p |−→ p | , (23) y − = 1 s [( √ s − p − p ) − |−→ p | − ( p ) − |−→ p | p ] , (24) y + = 1 s [( √ s − p − p ) − |−→ p | − ( p ) + 2 |−→ p | p ] , (25)where y is a dimensionless parameter defined as y = ( p − p − p ) /s with √ s = m t , and |−→ p | = q ( p ) − m B c . (26)The sum of the soft and hard sectors gives the total contribution of real corrections, i.e.,Γ Real = Γ
IRReal + Γ hardReal .With the real and virtual corrections, we then obtain the total decay width of t quarkto ¯ B c and Υ at the NLO accuracy of QCDΓ total = Γ Born + Γ
V irtual + Γ
Real + O ( α s ) . (27)In above expression, the decay width is UV and IR finite. In our calculation the FeynArts [16]was used to generate the Feynman diagrams, the amplitudes were generated by the FeynCalc[17], and the LoopTools [18] was employed to calculate the Passarino-Veltman integrations.The numerical integrations of the phase space were performed by the MATHEMATICA. IV. NUMERICAL RESULTS
To complete the numerical calculation, the following ordinarily accepted input parame-ters are taken into account: m b = 4 . , m c = 1 . , m t = 174 GeV , m W = 80 GeV , (28) ψ ¯ B c (0) = ψ ¯ B ∗ c (0) = R (0) √ π = 0 . / , (29) ψ LO Υ (0) = ψ LOη b (0) = R LO (0) √ π = 0 . / , (30) ψ NLO Υ (0) = ψ NLOη b (0) = R NLO (0) √ π = R LO (0) √ π − C F α s = 0 . / , (31)10 ABLE I: The decay widths of the processes t → ¯ B ∗ c + W + + c , t → ¯ B c + W + + c , t → Υ + W + + b and t → η b + W + + b at the tree level and with the NLO QCD corrections are presented in tworenormalization scale µ limits, those are 2 m c and m t for the first two processes and 2 m b and m t for the other two. t → ¯ B ∗ c + W + + c t → ¯ B c + W + + c t → Υ + W + + b t → η b + W + + bµ m c m t m c m t m b m t m b m t Γ LO NLO V tb = 1 . , G F = 1 . × − GeV − . (32)Here, V tb is the Cabibbo-Kobayashi-Maskawa(CKM) matrix element and G F is weak inter-action Fermi constant.In above numerical calculation inputs, the radial wave function at the origin for S-wave ¯ B ∗ c ( ¯ B c ) system is estimated by potential model [19], while the corresponding Υ( η b )nonperturbative parameter is determined from its electronic decay rate [8]. One loop resultof strong coupling constant is taken into account, i.e. α s ( µ ) = 4 π (11 − n f )Log( µ Λ QCD ) . (33)With the above preparation, one can readily obtain the decay widths of top quark to b ¯ c and b ¯ b mesons, as listed in Table I. To see the scale dependence of the LO and NLO results,the ratios Γ( µ ) / Γ(2 m c ) for b ¯ c system and Γ( µ ) / Γ(2 m b ) for b ¯ b system are showed in Figures6 and 7, respectively. Calculation tells that after including the NLO corrections, the energyscale dependence of the results is reduced, as expected. V. SUMMARY AND CONCLUSIONS
In this work we have calculated the decay widths of top quark to S-wave b ¯ c and b ¯ b boundstates at the NLO accuracy of perturbative QCD. Considering that there will be copious t ¯ t data in the near future at the LHC, our results are helpful to the study of the indirect11
20 40 60 80 100 120 140 160 1800.10.20.30.40.50.60.70.80.91.0 () / ( m c ) (GeV) NLO LO () / ( m c ) (GeV) NLO LO FIG. 6: The ratio Γ( µ ) / Γ(2 m c ) versus renormalization scale µ in t quark decays. The left diagramfor the b ¯ c spin-singlet state ¯ B c and the right diagram for the spin-triplet state ¯ B ∗ c .
20 40 60 80 100 120 140 160 1800.30.40.50.60.70.80.91.0 () / ( m b ) (GeV) NLO LO
20 40 60 80 100 120 140 160 1800.30.40.50.60.70.80.91.0 () / ( m b ) (GeV) NLO LO FIG. 7: The ratio Γ( µ ) / Γ(2 m b ) versus renormalization scale µ in t quark decays. The left diagramfor the b ¯ b spin-singlet state η b and the right diagram for the spin-triplet state Υ. production of these states. They may be also useful to the future study of NLO heavy quarkto b ¯ c and b ¯ b bound states fragmentation functions.Numerical results indicate that the NLO corrections greatly enhance the LO results for b ¯ b system, while slightly decrease the b ¯ c states production widths. The main reason forthis difference is that the NLO wave function for bottomonium is much larger than that ofLO one, while for the calculation of ¯ B c meson, the same wave function given by potentialmodel is used. Although from Table I, superficially the number of indirectly produced ¯ B c overshoots that of Υ, experimentally to detect the latter is much easier than the former.12ince top quark dominantly decays into b and W + final state with a width of 1.5 GeV orso, numerical results remind us that the Υ indirect production from top quark decay isdetectable, while it is hard to pin down the ¯ B c states by this way.The numerical calculation also shows that the next-to-leading order QCD correctionsto processes t → b ¯ c ( b ¯ b ) + W + + c ( b ) decrease the energy scale dependence of the decaywidths as expected, and hence the uncertainties in theoretical estimation. Future preciseexperiment on the concerned processes may provide a test on the theoretical framework forheavy quarkonium production and the reliability of perturbative calculation for them. Acknowledgments
This work was supported in part by the National Natural Science Foundation ofChina(NSFC) under the grants 10935012, 10928510, 10821063 and 10775179, by the CASKey Projects KJCX2-yw-N29 and H92A0200S2. [1] W. Hollik, in Proceedings of the XVI International Symposium on Lepton-Photon Interactions,Connell University, Ithaca, N.Y., Aug. 10-15 1993; M. Swartz, in Proceedings of the XVIInternational Symposium on Lepton-Photon Interactions, Connell University, Ithaca, N.Y.,Aug. 10-15 1993.[2] G. Altarelli, in Proceedings of Interational University School of Nuclear and Partical Physics:Substructures of Matter as Revealed with Electroweak Probes, Schladming, Ausria, 24 Feb -5 Mar. 1993.[3] G. Altarelli, CERN-TH-7319/94, talk at 1st International Conference on Phenomenology ofUnification: from Present to Future, Rome, Itali, 23-26 Mar 1994.[4] G. L. Kane, in Proceedings of the Workshop on High Energy Phenomenology, Mexico City,July 1-10, 1991.[5] F. Abe, et al. (CDF Collaboration), Phys. Rev. Lett. 74, 2626 (1995); S. Abachi, et al. (D0Collaboration), Phys. Rev. Lett. 74, 2632 (1995).[6] N. Kidonakis and R. Vogt, Int. J. Mod. Phys. A20, 3171, (2005); F. Hubaut, et al., ATLAScollaboration, hep-ex/0605029; V. Barger and R. J. Phillips, Preprint MAD/PH/789, 1993.
7] E. Braaten and T.C. Yuan, Phys. Rev. Lett. 71, 1673 (1993); ibid, Phys. Rev. D50, 3176(1994); C.-H. Chang and Y.-Q. Chen, Phys. Lett. B284, 127 (1992); ibid, Phys. Rev. D46,3845 (1992); Y.-Q. Chen, Phys. Rev. D48, 5181 (1993); T.C. Yuan, Phys. Rev. D50, 5664(1994).[8] E. Braaten, K. Cheung and T.C. Yuan, Phys. Rev. D48, 4230 (1993); ibid, Phys. Rev. D48,R5049 (1993).[9] K. Kolodzidj, A. Leike and R. Ruckl, Phys. Lett. B355, 337 (1995); C.-H. Chang, Y.-Q. Chen,Phys. Lett. B364, 78 (1995); C.-H. Chang, Y.-Q. Chen, Phys. Rev. D48, 4086 (1993); C.-H.Chang, Int. J. Mod. Phys. A21, 777 (2006); C.-H. Chang, J.-X. Wang and X.-G. Wu, Phys.Rev. D70, 114019 (2004); C.-H. Chang, C.-F. Qiao, J. X. Wang and X.-G. Wu, Phys. Rev.D71, 074012 (2005); C.-H. Chang, X.-G. Wu, Eur.Phys.J.C38, 267(2004); C.-H. Chang, C.Driouichi, P. Eerola and X.-G. Wu, Comput. Phys. Commun. 159, 192(2004); C.-H. Chang,J.-X. Wang and X.-G. Wu, Comput. Phys. Commun. 174, 241(2006).[10] I.I. Bigi, Yu L. Dokshitzer, V.A. Khoze, J.H. Kuhn and P. Zerwas, Phys. Lett. B181, 157(1986); L.H. Orr and J.L. Rosner, Phys. Lett. B246, 221 (1990).[11] C.-F. Qiao, C.-S. Li, K.-T. Chao, Phys. Rev. D54, 5606 (1996).[12] C.-H. Chang, J.-X. Wang, X.-G. Wu, Phys. Rev. D77, 014022 (2008).[13] X.-G. Wu, Phys. Lett. B671, 318 (2009).[14] G.T. Bodwin, E. Braaten, and G. P. Lepage, Phys. Rev. D51, 1125 (1995).[15] Y.-J. Zhang, Y.-J. Gao and K.-T. Chao, Phys. Rev. Lett. 96, 092001 (2006).[16] T. Hahn, Comput. Phys. Commun. 140, 418 (2001).[17] R. Mertig, M. B¨ohm, and A. Denner, Comput. Phys. Commun. 4, 345 (1991).[18] T. Hahn, M. Perez-Victoria Comput. Phys. Commun. 118, 153 (1999).[19] E.J. Eichten and C. Quigg, Phys. Rev. D52, 1726 (1995).7] E. Braaten and T.C. Yuan, Phys. Rev. Lett. 71, 1673 (1993); ibid, Phys. Rev. D50, 3176(1994); C.-H. Chang and Y.-Q. Chen, Phys. Lett. B284, 127 (1992); ibid, Phys. Rev. D46,3845 (1992); Y.-Q. Chen, Phys. Rev. D48, 5181 (1993); T.C. Yuan, Phys. Rev. D50, 5664(1994).[8] E. Braaten, K. Cheung and T.C. Yuan, Phys. Rev. D48, 4230 (1993); ibid, Phys. Rev. D48,R5049 (1993).[9] K. Kolodzidj, A. Leike and R. Ruckl, Phys. Lett. B355, 337 (1995); C.-H. Chang, Y.-Q. Chen,Phys. Lett. B364, 78 (1995); C.-H. Chang, Y.-Q. Chen, Phys. Rev. D48, 4086 (1993); C.-H.Chang, Int. J. Mod. Phys. A21, 777 (2006); C.-H. Chang, J.-X. Wang and X.-G. Wu, Phys.Rev. D70, 114019 (2004); C.-H. Chang, C.-F. Qiao, J. X. Wang and X.-G. Wu, Phys. Rev.D71, 074012 (2005); C.-H. Chang, X.-G. Wu, Eur.Phys.J.C38, 267(2004); C.-H. Chang, C.Driouichi, P. Eerola and X.-G. Wu, Comput. Phys. Commun. 159, 192(2004); C.-H. Chang,J.-X. Wang and X.-G. Wu, Comput. Phys. Commun. 174, 241(2006).[10] I.I. Bigi, Yu L. Dokshitzer, V.A. Khoze, J.H. Kuhn and P. Zerwas, Phys. Lett. B181, 157(1986); L.H. Orr and J.L. Rosner, Phys. Lett. B246, 221 (1990).[11] C.-F. Qiao, C.-S. Li, K.-T. Chao, Phys. Rev. D54, 5606 (1996).[12] C.-H. Chang, J.-X. Wang, X.-G. Wu, Phys. Rev. D77, 014022 (2008).[13] X.-G. Wu, Phys. Lett. B671, 318 (2009).[14] G.T. Bodwin, E. Braaten, and G. P. Lepage, Phys. Rev. D51, 1125 (1995).[15] Y.-J. Zhang, Y.-J. Gao and K.-T. Chao, Phys. Rev. Lett. 96, 092001 (2006).[16] T. Hahn, Comput. Phys. Commun. 140, 418 (2001).[17] R. Mertig, M. B¨ohm, and A. Denner, Comput. Phys. Commun. 4, 345 (1991).[18] T. Hahn, M. Perez-Victoria Comput. Phys. Commun. 118, 153 (1999).[19] E.J. Eichten and C. Quigg, Phys. Rev. D52, 1726 (1995).