NLO QCD corrections to B c -pair production in photon-photon collision
aa r X i v : . [ h e p - ph ] J u l NLO QCD corrections to B c -pair production inphoton-photon collision Zi-Qiang Chen † , Hao Yang ‡ and Cong-Feng Qiao , ∗ School of Physics, University of Chinese Academy of Science,Yuquan Road 19A, Beijing 10049 CAS Key Laboratory of Vacuum Physics, Beijing 100049, China
Abstract
The B c meson pair, including pairs of pseudoscalar states and vector states, productions inhigh energy photon-photon interaction are investigated at the next-to-leading order (NLO)accuracy in the nonrelativistic quantum chromodynamics (NRQCD) factorization formalism.The corresponding cross sections at the future e + e − colliders with √ s = 250 GeV and 500GeV are evaluated. Numerical result indicates that the inclusion of the NLO corrections shallgreatly suppress the scale dependence and enhance the prediction reliability. In addition tothe phenomenological meaning, the NLO QCD calculation of this process subjects to certaintechnical issues, which are elucidated in details and might be applicable to other relevantinvestigations.PACS numbers: 12.38.Bx, 12.39.Jh, 14.40.Pq, 14.70.Bh † [email protected] ‡ [email protected] ∗ [email protected], corresponding author . INTRODUCTION As the only heavy meson consisting of two heavy quarks with different flavors, the B c meson is of great interest in both experiment and theory. Study of its productionand decays may enrich our knowledge on the properties of double heavy meson andthe nature of perturbative QCD (pQCD). The ground state of B c meson, B + c (1 S ), wasdiscovered by CDF Collaboration [1, 2] in 1998. And its excited state B + c (2 S ) wasobserved by ATLAS [3] and CMS [4] Collaborations in 2014 and 2019 respectively.Due to the large masses of bottom and charm quarks, the production of heavyquark pair can be described by pQCD, while the hadronization process can be factoredby using the NRQCD factorization formalism [5]. For inclusive B c meson production,various investigations have been carried out, including the direct production through pp [6–9], e + e − [10, 11], γγ [12, 13] and ep [14, 15] collisions, and the indirect productionthrough top quark [16, 17], Z boson [18–21], W boson [22–24] and Higgs boson [25]decays.Within QCD and quantum electromagnetic dynamics (QED), the B + c meson is pro-duced in accompany with an additional b ¯ c pair, which is also possible to form another b ¯ c meson, namely the B c -pair may exclusively produced. Generally speaking, the ex-periment measurement of exclusive process possesses a relative high precision, whichis required in exploring the properties of QCD and hadrons. In the literature, various B c -pair production processes have been investigated, including in pp [26, 27], e + e − [28–30] and γγ [26] collisions. We notice that in Ref.[26] the leading order (LO) analysison B c -pair production in photon-photon collision was performed, however with only B + c + B − c (pseudoscalar-pseudoscalar, PP) and B ∗ + c + B ∗− c (vector-vector, VV) config-urations being considered. In this work, for the sake of completeness we first repeatthe LO calculation in [26] and then calculate the LO B c -pair production in B + c + B ∗− c (pseudoscalar-vector, PV) and B ∗ + c + B − c (vector-pseudoscalar, VP) configurations . In The PV and VP production are related by a charge-conjugation transformation. Their cross sectionare exactly the same. B c represents for both pseudoscalar B c and vector B ∗ c , the lattermay overwhelmingly decay to the pseudoscalar state, unless specifically mentioned.The rest of the paper is organized as follows. In section II we present the primaryformulae employed in the calculation. In section III, some technical details in theanalytical calculation are given. In section IV, the numerical evaluation for concernedprocesses is performed. The last section is remained for summary. II. FORMULATION
According to NRQCD factorization formalism, the cross section of B c -pair produc-tion via photon-photon fusion can be formulated as d ˆ σ ( γ + γ → B + c + B − c ) = | ψ (0) | s X |M ( γ + γ → [ c ¯ b ] + [ b ¯ c ]) | d PS , (1)where ψ (0) is the wave function of B c meson at the origin, ˆ s is the center-of-massenergy square of two colliding photons, P sums over the polarizations and colors ofthe initial and final particles, comes from the spin average of the initial γγ states, M ( γ + γ → [ c ¯ b ] + [ b ¯ c ]) is the corresponding partonic amplitude, d PS stands for thetwo-body phase space.The partonic amplitude can be computed by using the covariant projection operatormethod. At the leading order of the relative velocity expansion, it is legitimate to take m B c = m b + m c , p B c = p c + p b = (1 + m c m b ) p b . The spin and color projection operatorhas the form Π( n ) = 12 √ m B c ǫ ( n )( /p B c + m B c ) ⊗ (cid:18) c √ N c (cid:19) , (2)where ǫ ( S ) = γ , ǫ ( S ) = /ǫ , and ǫ represents the polarization vector of B ∗ c meson.The 1 c stands for the unit color matrix, and N c = 3 for the number of colors in QCD.The photon-photon scattering may be achieved in high energy e + e − collider likethe Large Electron-Positron Collider (LEP), the Circular Electron Positron Collider(CEPC) and the International Linear Collider (ILC), or even in hadron collider like3he Large Hadron Collider (LHC). Here we focus only on the e + e − collision case,where the initial photon can be generated by the bremsstrahlung or by the laser backscattering (LBS) effect. The cross section are then formulated as dσ ( e + + e − → e + + e − + B + c + B − c ) = Z dx dx f γ ( x ) f γ ( x ) dσ ( γ + γ → B + c + B − c ) , (3)where f γ ( x ) is the photon distribution with fraction x of the beam energy.Imposing transverse momentum cut p − T < p T < p + T and rapidity cut | y | < y c oneach B c meson, the formula for total cross section is then σ ( e + + e − → e + + e − + B + c + B − c )= 1256 π ( θ (cid:0) − ln m − T √ s (cid:1) Z min { , ln m + T √ s } ln m − T √ s dX Z min { y + T ,y c } max {− y + T , − y c } dy ∗ sech y ∗ E X |M| Z min { y c − y ∗ , − X } max {− y c + y ∗ ,X } dy x f γ ( x ) x f γ ( x ) + θ (cid:0) − ln m + T √ s (cid:1) Z min { , ln( m + T √ s cosh y c ) } ln m + T √ s dX (cid:18) Z min { y + T ,y c } y − T dy ∗ + Z − y − T max {− y + T , − y c } dy ∗ (cid:19) sech y ∗ E X |M| Z min { y c − y ∗ , − X } max {− y c + y ∗ ,X } dy x f γ ( x ) x f γ ( x ) ) , (4)with X = 12 ln( x x ) , y = 12 ln x x ,m ± T = q m B c + p ± T ,y ± T = 12 ln E + q E − m ∓ T E − q E − m ∓ T . (5)Here, √ s is the collision energy of e + e − collider, E = √ sx x / y ∗ = y − y arerespectively the energy and rapidity of B c meson in the photon-photon center-of-masssystem, θ ( x ) means the unit step function.The spectrum of bremsstrahlung photon is well formulated in the Weizsacker-Williams approximation (WWA) as [31] f γ ( x ) = α π (cid:20) − x ) x log (cid:18) Q Q (cid:19) + 2 m e x (cid:18) Q − Q (cid:19)(cid:21) , (6)4 .0 0.2 0.4 0.6 0.8 1.0 x012345 f γ ( x ) LBSWWA s =
250 GeV
FIG. 1: The spectra of WWA photon and LBS photon at √ s = 250 GeV. where Q = m e x / (1 − x ) and Q = Q + ( θ c √ s/ (1 − x ) with x = E γ /E e , θ c is the experimental angular cut which taken to be 32 mrad here. For the LBS photon,the spectrum is expressed as [32] f γ ( x ) = 1 N (cid:20) − x + 11 − x − r (1 − r ) (cid:21) , (7)where r = xx m (1 − x ) and the normalization factor N = (cid:18) − x m − x m (cid:19) log(1 + x m ) + 12 + 8 x m − x m ) . (8)Here x m ≃ .
83 [33] and the energy fraction x of photon is restricted in 0 ≤ x ≤ x m / (1 + x m ). The behaviors of WWA photon and LBS photon are quite different, theirspectra at √ s = 250 GeV are shown in Fig.1. III. ANALYTICAL CALCULATION
The typical tree-level and one-loop Feynman diagrams for the partonic processes areshown in Fig.2. The momenta and the polarization vectors of incoming and outgoingparticles are denoted as: γ ( p , ǫ ) + γ ( p , ǫ ) → [ c ¯ b ]( k , ǫ ) + [¯ cb ]( k , ǫ ) . (9)5 b cb c b ( c ) ( d ) ( e )( f ) cb ( g ) c b ( h ) cb ( i ) u, d, s, c, b ( j ) c b ( k ) cb ( m ) cb ( o )( n ) cbc b ( a ) ( b )( l ) c bc b cb FIG. 2: Typical Feynman diagrams of the partonic processes for B c -pair production. (a) and(b): tree-level diagrams; (c) and (d): self-energy corrections; (e) and (f): vertex corrections;(g)-(i): box diagrams; (j) and (k): pentagon diagrams; (l) and (m): hexagon diagrams; (n)and (o): diagrams of counter terms. Here, initial and final state particles are all on their mass shells: p = p = 0 and k = k = m B c . The polarization vectors satisfy the constraints: ǫ · ǫ ∗ = ǫ · ǫ ∗ = ǫ · ǫ ∗ = ǫ · ǫ ∗ = − p · ǫ = p · ǫ = k · ǫ = k · ǫ = 0.To proceed the calculation, we notice working in the photon-photon center-of-mass system is convenient. By introducing the orthonormal four-vector base: n =(1 , , , n = (0 , , , n = (0 , , ,
0) and n = (0 , , , p = E ( n + n ) , p = E ( n − n ) ,k = E ( n + r y n + r z n ) , k = E ( n − r y n − r z n ) . (10)and ǫ (1)1 = n , ǫ (2)1 = n , ǫ (1)2 = n , ǫ (2)2 = n ;6 (1)3 = n , ǫ (2)3 = r z n − r y n p r y + r z , ǫ (3)3 = ( r y + r z ) n + r y n + r z n r m p r y + r z ; ǫ (1)4 = n , ǫ (2)4 = r z n − r y n p r y + r z , ǫ (3)4 = ( r y + r z ) n − r y n − r z n r m p r y + r z . (11)Here, E = √ sx x / r y = k y /E , r z = k z /E , r m = m B c /E , and the on-shellcondition constrains r y + r z + r m = 1. Then, the helicity amplitudes can be readilycalculated through M ij PP = A µν PP ǫ ( i )1 µ ǫ ( j )2 ν , M ijk PV = A µνρ PV ǫ ( i )1 µ ǫ ( j )2 ν ǫ ( k )3 ρ , M ijk VP = A µνρ VP ǫ ( i )1 µ ǫ ( j )2 ν ǫ ( k )4 ρ , M ijkl VV = A µνρσ VV ǫ ( i )1 µ ǫ ( j )2 ν ǫ ( k )3 ρ ǫ ( l )4 σ . (12)The tree-level calculation is straightforward, however the full analytic expressions ofhelicity amplitudes are still too lengthy to present in the mainbody of text. Consideringof the symmetric property in amplitudes, we present the LO results in Appendix.In the computation of one-loop amplitudes, the conventional dimensional regular-ization with D = 4 − ǫ is adopted to regularize the ultraviolet (UV) and infrared (IR)singularities. The IR singularities are canceled each other and the UV singularitiesare removed by renormalization procedure. The renormalization constants include Z , Z m , Z and Z g , corresponding to heavy quark field, heavy quark mass, gluon fieldand strong coupling constant, respectively. We define Z and Z m in the on-shell (OS)scheme, Z and Z g in the modified minimal-subtraction (MS) scheme. The correspond-ing counterterms are δZ OS2 = − C F α s π (cid:20) ǫ UV + 2 ǫ IR − γ E + 3 ln 4 πµ m + 4 (cid:21) ,δZ OS m = − C F α s π (cid:20) ǫ UV − γ E + ln 4 πµ m + 43 (cid:21) ,δZ MS3 = α s π ( β − C A ) (cid:20) ǫ UV − γ E + ln(4 π ) (cid:21) ,δZ MS g = − β α s π (cid:20) ǫ UV − γ E + ln(4 π ) (cid:21) . (13)7ere, µ is the renormalization scale, γ E is the Euler’s constant; m stands for m c and m b accordingly; β = (11 / C A − (4 / T f n f is the one-loop coefficient of QCD betafunction, n f is the number of active quarks which taken to be 5 in our calculation; C A = 3, C F = 4 / T F = 1 / δZ terms cancel with each other.For reference, we provide the analytic results for one-loop amplitudes as supplemen-tary files attached to the arXiv preprint. In the NLO calculation, the Mathematicapackage FeynArts [34] is used to generate Feynman diagrams and Feynman amplitudes;FeynCalc [35, 36] and FORM [37, 38] are used to perform algebraic calculation; Thepackage FIRE [39, 40] is employed to reduce the Feynman integrals to typical masterintegrals A , B , C and D , which are numerically evaluated by LoopTools [41]. IV. NUMERICAL RESULTS
In numerical analysis, the formula (4) is employed with |M| ≃ |M tree | for theLO calculation and |M| ≃ |M tree | + 2Re( M loop M ∗ tree ) for the NLO calculation. Therapidity and p T cuts, | y | < < p T <
40 GeV, are imposed on each B c meson.Other inputs in numerical evaluation go as follows: α = 1 / . , m e = 0 .
511 Mev , m c = 1 . ,m b = 4 . , | ψ (0) | = 0 .
174 GeV . (14)Here, the B c wave function at the origin is estimated from the S − S splitting [42] | ψ (0) | = 9 m b m c πα s ( M B ∗ c − M B c ) (15)with the lattice calculation result on M B ∗ c − M B c = 53 MeV [43].The two-loop strong coupling of α s ( µ )4 π = 1 β L − β ln Lβ L (16)is employed in the NLO calculation, in which, L = ln( µ / Λ ), β = (34 / C A − C F T F n f − (20 / C A T F n f , with n f = 5 and Λ QCD = 210 MeV adopted here [44].8ote, for LO calculation, the one-loop formula of the running coupling constant isused.Considering in future the e + e − collider like CEPC and ILC might run at center-of-mass energies √ s = 250 GeV and √ s = 500 GeV respectively, we numerically evaluatethe B c -pair production via WWA and LBS schemes at these two energies. Taking thesame inputs, we can numerically repeat the LO double pseudoscalar B c productionresult in [26]. The full NLO results are presented in Fig.3, Fig.4 and Fig.5. Note,because the cross sections for PV production and VP production are exactly the same,only the PV production results are illustrated.The total cross sections versus r are shown in Fig.3 with µ = r q m B c + p T . Weobserve that, in comparison with the LO contribution, the LO plus NLO cross sectionsare suppressed, as are their dependences on the renormalization scale. As the √ s rises from 250 GeV to 500 GeV, the B c -pair production rates increase in the WWAmechanism while decrease in the LBS mechanism. This may be understood from thedifferent behaviors of WWA and LBS mechanisms, as shown in Fig.1. The WWAphotons are more likely to be produced with small momentum fraction x , while theLBS photons tend to be more energetic. Moreover, the partonic cross section ˆ σ ( γ + γ → B + c + B − c ) decreases with the increase of incident photons’ center-of-mass energy.The differential cross sections as functions of p T , the transverse momentum of oneof the two B c mesons, are shown in Fig.4. It can be seen that as √ s increases from250 GeV to 500 GeV, the yields of WWA processes increase slightly, while the yieldsof LBS processes decrease, evidently in small p T region. Since the LBS photons aregenerally more energetic than the WWA photons, the produced B c pairs tend to havelarger transverse momenta, which shall lead to a flatter p T distribution.The differential cross sections as functions of ∆ y , the rapidity difference betweentwo produced B c mesons, are shown in Fig.5. Note, due to | ∆ y | = 2 | y ∗ | , the | ∆ y | distribution is equivalent to the | y ∗ | distribution, where y ∗ is the rapidity of B c mesonin the photon-photon center-of-mass frame. For PP and VV production, the B c pairsare more likely to be produced around the y ∗ = 0 region, while for the PV or VP9 .6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 r0.150.200.250.300.350.40 σ PP ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
250 GeV, PP (a) σ PP ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
500 GeV, PP (b) σ PV ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
250 GeV, PV (c) σ PV ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
500 GeV, PV (d) σ VV ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
250 GeV, VV (e) σ VV ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
500 GeV, VV (f)
FIG. 3: The LO and NLO cross sections versus r , where r = µ q m Bc + p T . (a) PP, 250 GeV;(b) PP, 500 GeV; (c) PV, 250 GeV; (d) PV, 500 GeV; (e) VV, 250 GeV; and (f) VV, 500GeV. production, the peak is located close to y ∗ = 0 .
6. Since the large energy may leadto large y ∗ , for the same reason as explained in p T distribution, the LBS productiondistributions are flatter than the WWA ones.10
10 15 20 25 30 35 40 p T ( GeV ) - - - d σ PP / dp T ( fb / GeV ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
250 GeV, PP (a) p T ( GeV ) - - σ PP / dp T ( fb / GeV ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
500 GeV, PP (b) p T ( GeV ) - - - - σ PV / dp T ( fb / GeV ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
250 GeV, PV (c) p T ( GeV ) - - - σ PV / dp T ( fb / GeV ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
500 GeV, PV (d) p T ( GeV ) - - σ VV / dp T ( fb / GeV ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
250 GeV, VV (e) p T ( GeV ) - - σ VV / dp T ( fb / GeV ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
500 GeV, VV (f)
FIG. 4: The LO and NLO differential cross sections versus p T , the transverse momentum ofone of the two B c mesons. The renormalization scale µ = q m B c + p T . (a) PP, 250 GeV; (b)PP, 500 GeV; (c) PV, 250 GeV; (d) PV, 500 GeV; (e) VV, 250 GeV; and (f) VV, 500 GeV. In Ref.[30], the production of B c pairs in e + e − annihilation via virtual γ ∗ and Z ∗ are investigated at the NLO QCD accuracy. At large √ s , say √ s >
160 GeV, the crosssections are less than 10 − fb, which are four to six orders of magnitude smaller thanthe cross sections of the processes considered here. It means that at high energy e + e − collider, photon-photon collision turns out to be the dominant mechanism for B c -pair11 | Δ y | - σ PP / d | Δ y | ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
250 GeV, PP (a) | Δ y | - σ PP / d | Δ y | ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
500 GeV, PP (b) | Δ y | - σ PV / d | Δ y | ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
250 GeV, PV (c) | Δ y | × - × - σ PV / d | Δ y | ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
500 GeV, PV (d) | Δ y | σ VV / d | Δ y | ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
250 GeV, VV (e) | Δ y | σ VV / d | Δ y | ( fb ) LBS, NLOLBS, LOWWA, NLOWWA, LO | y |<
2, 2 < p T <
40 GeV s =
500 GeV, VV (f)
FIG. 5: The LO and NLO differential cross sections versus ∆ y , the rapidity difference betweenthe two B c mesons. The renormalization scale µ = q m B c + p T . (a) PP, 250 GeV; (b) PP,500 GeV; (c) PV, 250 GeV; (d) PV, 500 GeV; (e) VV, 250 GeV; and (f) VV, 500 GeV. production.Moreover, about the numerical results, there are some points remarkable:1) At √ s =250 GeV, as can be seen from Fig.1, the WWA and LBS spectra meet at12bout x ∼ .
07, and hence the WWA and LBS cross sections tend to be comparable inthe corresponding kinematic region (i.e. in about 2 < p T < < | ∆ y | < p T andsmall | ∆ y | regions dominate the production, the WWA and LBS cross sections arehence almost identical, as shown in Fig.3(e).2) To investigate the convergence of perturbative expansion, we define a measure R = | σ NLO − σ LO σ LO | . For PP, PV and VV productions, we find0 < R PP < . , . < R PV < . , . < R VV < . , (17)which are compatible with the results in Ref.[30], where the B c -pair production via e + e − annihilation was studied. Furthermore, the NLO effect is more significant in thesmall p T region, suggests a resummation for the logarithms of p T ˆ s , which is beyond thescope of this work.3) In the numerical calculation, the strong coupling constant is parameterized by Λ QCD as in Eq.(16). Whereas, it is noteworthy that nowadays an alternative approach iswidely accepted, in which the initial value of α s at a well experimentally measuredpoint is adopted, usually at M Z , rather than Λ QCD . Then, for fixed n f , one may useevolution equation [45]4 πα s ( µ ) − β β ln (cid:18) πα s ( µ ) + β β (cid:19) = 4 πα s ( M Z ) − β β ln (cid:18) πα s ( M Z ) + β β (cid:19) + β ln µ M Z (18)to run the α s to where interested in. We evaluate the cross sections as well employingthis parameterizing scheme with input α s ( M Z ) = 0 . QCD determined by α s ( M Z ), i.e. Λ LOQCD = 88 MeV, Λ
NLOQCD = 228MeV. In all, the discrepancy in strong coupling constant parameterization between twoschemes may be somehow remedied by adjusting the value of Λ
QCD , but the evolutionscheme is recommended. 13 . SUMMARY
In this work we investigated the B c -pair production in high energy photon-photonfusion at the NLO accuracy in the NRQCD factorization framework. Various of S -wave B c states, including configurations of PP, PV, VP and VV, were taken intoaccount. Considering the leading order results for B c -pair production in PV and VPconfigurations are still missing in the literature, we calculated them and provided theanalytic results. The total cross section as well as p T and ∆ y distributions in e + e − collider with √ s =250 GeV and √ s = 500 GeV were evaluated and presented in figures.The numerical results show that with the NLO corrections, the LO cross sectionsare suppressed, and their dependence on renormalization scale are reduced evidently.By comparing with the results in Ref.[30], where the B c -pair production in e + e − anni-hilation via virtual γ ∗ and Z ∗ were investigated, we may conclude that at large e + e − collision energy, say √ s >
160 GeV, photon-photon collision will be the dominantsource of B c -pair production.The NLO calculation of the concerned processes is somewhat time consuming andcomputer resource exhausting. To fulfill this work, a ”divide-and-conquer” strategywas employed. Instead of squaring the amplitude and summing over spins, we calcu-lated the helicity amplitudes separately, which makes this tedious calculation workable.Moreover, it shows that the symmetries remain in the helicity amplitudes may greatlyreduce the number of independent amplitudes, as illustrated in Appendix.Last, the concerned processes involve a number of momenta and polarization vectorsof the external particles, by introducing auxiliary vectors, the base n , n , n and n , thenumber of independent Lorentz vectors reduces to 4, which facilitate the computationof Feynman integrals. We think the technical strategy employed in this work might beapplicable to the studies of some other relevant processes.Note added: when this work was finished and the manuscript was finalizing, thereappeared a study on the web about the B c -pair production in photon-photon collisonwith the relativistic corrections [46]. 14 cknowledgments This work was supported in part by the National Natural Science Foundation ofChina(NSFC) under the Grants 11975236 and 11635009. [1] F. Abe et al. [CDF], Phys. Rev. Lett. , 2432-2437 (1998) [arXiv:hep-ex/9805034 [hep-ex]].[2] F. Abe et al. [CDF], Phys. Rev. D , 112004 (1998) [arXiv:hep-ex/9804014 [hep-ex]].[3] G. Aad et al. [ATLAS], Phys. Rev. Lett. , no.21, 212004 (2014) [arXiv:1407.1032[hep-ex]].[4] A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. , no.13, 132001 (2019)[arXiv:1902.00571 [hep-ex]].[5] G. T. Bodwin, E. Braaten and G. Lepage, Phys. Rev. D , 1125-1171 (1995); ErratumPhys. Rev. D , 5853 (1997) [arXiv:hep-ph/9407339 [hep-ph]].[6] C. H. Chang and Y. Q. Chen, Phys. Rev. D , 4086-4091 (1993)[7] C. H. Chang, Y. Q. Chen, G. P. Han and H. T. Jiang, Phys. Lett. B , 78-86 (1995)[arXiv:hep-ph/9408242 [hep-ph]].[8] A. Berezhnoy, A. Likhoded and O. Yushchenko, Phys. Atom. Nucl. , 709-713 (1996)[arXiv:hep-ph/9504302 [hep-ph]].[9] K. Kolodziej, A. Leike and R. Ruckl, Phys. Lett. B , 337-344 (1995)[arXiv:hep-ph/9505298 [hep-ph]].[10] Z. Yang, X. G. Wu and X. Y. Wang, Comput. Phys. Commun. , 2848-2855 (2013)[arXiv:1305.4828 [hep-ph]].[11] X. C. Zheng, C. H. Chang, T. F. Feng and Z. Pan, Sci. China Phys. Mech. Astron. ,no.3, 031012 (2018) [arXiv:1701.04561 [hep-ph]].[12] A. Berezhnoy, A. Likhoded and M. Shevlyagin, Phys. Lett. B , 351-355 (1995)[arXiv:hep-ph/9408287 [hep-ph]].[13] K. Kolodziej, A. Leike and R. Ruckl, Phys. Lett. B , 219-225 (1995) arXiv:hep-ph/9412249 [hep-ph]].[14] A. Berezhnoy, V. Kiselev and A. Likhoded, Phys. Atom. Nucl. , 252-259 (1998)[arXiv:hep-ph/9710429 [hep-ph]].[15] H. Y. Bi, R. Y. Zhang, H. Y. Han, Y. Jiang and X. G. Wu, Phys. Rev. D , no.3,034019 (2017) [arXiv:1612.07990 [hep-ph]].[16] C. F. Qiao, C. S. Li and K. T. Chao, Phys. Rev. D , 5606-5610 (1996)[arXiv:hep-ph/9603275 [hep-ph]].[17] P. Sun, L. P. Sun and C. F. Qiao, Phys. Rev. D , 114035 (2010) [arXiv:1003.5360[hep-ph]].[18] C. H. Chang and Y. Q. Chen, Phys. Rev. D , 3845 (1992); Erratum Phys. Rev. D ,6013 (1994).[19] V. Kiselev, A. Likhoded and M. Shevlyagin, Phys. Atom. Nucl. , 689-699 (1994)[arXiv:hep-ph/9401348 [hep-ph]].[20] C. F. Qiao, L. P. Sun and R. L. Zhu, JHEP , 131 (2011) [arXiv:1104.5587 [hep-ph]].[21] J. Jiang, L. B. Chen and C. F. Qiao, Phys. Rev. D , no.3, 034033 (2015)[arXiv:1501.00338 [hep-ph]].[22] C. F. Qiao, L. P. Sun, D. S. Yang and R. L. Zhu, Eur. Phys. J. C , 1766 (2011)[arXiv:1103.1106 [hep-ph]].[23] X. C. Zheng, C. H. Chang, X. G. Wu, J. Zeng and X. D. Huang, Phys. Rev. D ,no.3, 034029 (2020) [arXiv:1911.12531 [hep-ph]].[24] Z. Q. Chen, H. Yang and C. F. Qiao, Phys. Rev. D , no.3, 036009 (2020)[arXiv:1912.02140 [hep-ph]].[25] J. Jiang and C. F. Qiao, Phys. Rev. D , no.5, 054031 (2016) [arXiv:1512.01327 [hep-ph]].[26] S. Baranov, Phys. Rev. D , 2756-2759 (1997)[27] R. Li, Y. J. Zhang and K. T. Chao, Phys. Rev. D , 014020 (2009) [arXiv:0903.2250[hep-ph]].[28] V. Kiselev, Int. J. Mod. Phys. A , 465-476 (1995)
29] A. Karyasov, A. Martynenko and F. Martynenko, Nucl. Phys. B , 36-51 (2016)[arXiv:1604.07633 [hep-ph]].[30] A. Berezhnoy, A. Likhoded, A. Onishchenko and S. Poslavsky, Nucl. Phys. B , 224-242 (2017) [arXiv:1610.00354 [hep-ph]].[31] S. Frixione, M. L. Mangano, P. Nason and G. Ridolfi, Phys. Lett. B , 339-345 (1993)[arXiv:hep-ph/9310350 [hep-ph]].[32] I. Ginzburg, G. Kotkin, V. Serbo and V. I. Telnov, Nucl. Instrum. Meth. , 47-68(1983)[33] V. I. Telnov, Nucl. Instrum. Meth. A , 72-92 (1990)[34] T. Hahn, Comput. Phys. Commun. , 418-431 (2001) [arXiv:hep-ph/0012260 [hep-ph]].[35] R. Mertig, M. Bohm and A. Denner, Comput. Phys. Commun. , 345-359 (1991)[36] V. Shtabovenko, R. Mertig and F. Orellana, Comput. Phys. Commun. , 432-444(2016) [arXiv:1601.01167 [hep-ph]].[37] J. Vermaseren, Nucl. Phys. B Proc. Suppl. , 19-24 (2008) [arXiv:0806.4080 [hep-ph]].[38] J. Kuipers, T. Ueda, J. Vermaseren and J. Vollinga, Comput. Phys. Commun. ,1453-1467 (2013) [arXiv:1203.6543 [cs.SC]].[39] A. Smirnov, JHEP , 107 (2008) [arXiv:0807.3243 [hep-ph]].[40] A. V. Smirnov, Comput. Phys. Commun. , 182-191 (2015) [arXiv:1408.2372 [hep-ph]].[41] T. Hahn and M. Perez-Victoria, Comput. Phys. Commun. , 153-165 (1999)[arXiv:hep-ph/9807565 [hep-ph]].[42] E. J. Eichten and C. Quigg, Phys. Rev. D , 5845-5856 (1994) [arXiv:hep-ph/9402210[hep-ph]].[43] E. B. Gregory, C. T. H. Davies, I. D. Kendall, J. Koponen, K. Wong, E. Follana,E. Gamiz, G. P. Lepage, E. H. Muller, H. Na and J. Shigemitsu, Phys. Rev. D ,014506 (2011) [arXiv:1010.3848 [hep-lat]].[44] M. Tanabashi et al. [Particle Data Group], Phys. Rev. D , no.3, 030001 (2018)
45] A. Deur, S. J. Brodsky and G. F. de Teramond, Prog. Part. Nucl. Phys. , 1-74 (2016)[arXiv:1604.08082 [hep-ph]].[46] A. E. Dorokhov, R. N. Faustov, A. P. Martynenko and F. A. Martynenko,[arXiv:2005.06053 [hep-ph]]. Appendix
For PP, PV (or VP) and VV production, there are 4, 12 and 36 helicity amplitudesrespectively, whereas, half of them are zero. The nonzero helicity amplitudes are M , M , M , VP , M , VP , M , VP , M , VP , M , VP , M , VP ; M , M , M , M , M , M , M , M , M , M , M , M , M , M , M , M , M , M . (19)The processes of PV and VP productions are correlated in charge-conjugation trans-formation, their cross sections should be exactly the same. According convention (10)and (11), the amplitudes satisfy M ijk PV = ±M ijk VP , (20)where the plus sign corresponds to { i, j, k } = { , , } and { , , } , the minus signcorresponds to other cases. In addition, the helicity amplitudes satisfy also M = M (cid:12)(cid:12) k z →− k z , M = −M (cid:12)(cid:12) k z →− k z ; M = M (cid:12)(cid:12) k z →− k z , M = −M (cid:12)(cid:12) k z →− k z ; M = −M , M = −M , M = −M (cid:12)(cid:12) k z →− k z = −M (cid:12)(cid:12) k z →− k z = M , M = −M (cid:12)(cid:12) k z →− k z = M (cid:12)(cid:12) k z →− k z = −M . (21)The analytical expressions for helicity amplitudes can be classified in photon-quark18oupling, as M = 8 C A C F m B c π αα s E m b m c " e c f − e c e b ( f + f ) + e b f + X i = u,d,s e i f , (22)where e q represents the electric charge number of quark q , i.e. e c = e u = , e d = e s = e b = − . The coefficients f and f , f and f are related as per m c ↔ m b exchange: f PP , (cid:12)(cid:12) m c ↔ m b = f PP , , f PP , (cid:12)(cid:12) m c ↔ m b = f PP , ; f PV , (cid:12)(cid:12) m c ↔ m b = − f PV , , f PV , (cid:12)(cid:12) m c ↔ m b = − f PV , ; f VP , (cid:12)(cid:12) m c ↔ m b = − f VP , , f VP , (cid:12)(cid:12) m c ↔ m b = − f VP , ; f VV , (cid:12)(cid:12) m c ↔ m b = f VV , , f VV , (cid:12)(cid:12) m c ↔ m b = f VV , . (23)For the tree amplitudes, f is zero. The analytical results for other coefficients are f , = − r (1 − r z ) r − − r y ( r − − r z ) + r ( rr y − r +3) r − − r y (1 − r z ) ,f , = 1 − r y (1 − r z ) ,f , = − r y (2 r r y − r +4 r − − r z ) − r (1 − r z ) r − + r y (4 r − r r y − r +1)( r − − r z ) + r (3 rr y − r +3) r − ,f , = r − r − r y − r z − r y (2 r r y − r − rr y +4 r − − r z ) + 1; f , = ir m r y r z (cid:0) r ( r − − r z ) + − r z ) (cid:1) ,f , = − i r m r y r z (1 − r z ) ,f , = ir m r y √ r y + r z (cid:0) r + rr y − r +2( r − − r z ) + rr y + rr z − r +1)(1 − r z ) − rr − (cid:1) ,f , = ir m r y √ r y + r z (cid:0) rr y + rr z − r − r y − r z )(1 − r z ) + r − r z (cid:1) ,f , = i √ r y + r z (cid:0) − rr y (2 r − r y − r − − r z ) + r y ( rr y +2 rr z − r z )(1 − r z ) − r (2 r y +1) r − + r (1 − r z ) r − (cid:1) ,f , = i √ r y + r z (cid:0) rr y +2 rr z − r y − r z )(1 − r z ) − r − − r z (cid:1) ,f , = ir m r y r z (cid:0) r − − r z ) − r ( r − − r z ) (cid:1) ,f , = i r m r y r z − r (1 − r z ) ; f , = r (1 − r z ) r − − r − r y − r − − r z ) − r ( rr y +1) r − + r y (1 − r z ) ,f , = − r y (1 − r z ) + − r z − , , = r y + r z (cid:0) − r r y + r r y + rr y +3 r + r y − r − − r (1 − r z ) r − + r (1 − r z )(2 rr y + r +1) r − + ( r y +1)(2 r + r y − r − − r z ) + r y − r y +1) r y (1 − r z ) (cid:1) ,f , = r y + r z (cid:0) − r y (1 − r z ) + r y +1)1 − r z − r y − r z − (cid:1) ,f , = r m r y r z r y + r z (cid:0) r ( r y +1)( r − − r z ) − rr − + r y − r y +1)(1 − r z ) (cid:1) ,f , = r m r y r z r y + r z (cid:0) − r z − − r z ) (cid:1) ,f , = r y + r z (cid:0) − r r y +3 r r y +2 r − rr y − r + r y r − − r (1 − r z ) r − + ( r y +1) r y ( r − − r z ) + r (1 − r z )(2 rr y +3 r − r − + r y − r y +1) r y (1 − r z ) (cid:1) ,f , = r y + r z (cid:0) r y − r z − r y (1 − r z ) − r y − r z (cid:1) ,f , = r m √ r y + r z (cid:0) − r ( r y +1)( r − − r z ) − rr y + r z )(1 − r z ) + rr − (cid:1) ,f , = r m √ r y + r z (cid:0) − rr y − r y − − r z ) − − r z (cid:1) ,f , = r m r y √ r y + r z (cid:0) r (2 r − r y − r − − r z ) − rr y + rr z + r − − r z ) + rr − (cid:1) ,f , = r m r y √ r y + r z (cid:0) r − − r z − rr y + rr z + r − r y − − r z ) (cid:1) ,f , = r y (2 r r y − − r z ) + r (1 − r z ) r − − r r y − r r y − r r y − rr y − r + r y +2( r − − r z ) − r (3 rr y +1) r − ,f , = r y (2 r r y − rr y − − r z ) − r r y − rr y − r y − − r z − ,f , = r y + r z (cid:0) r y (2 r r y +2 r r y − r y +1)(1 − r z ) − r (1 − r z ) r − + r r y − r r y − r r y − rr y + r + r y − r − − r r y +4 r r y − r r y − r r y − r r y − rr y +2 r + r y − r y − r − − r z ) + r (1 − r z )(4 rr y + r +1) r − (cid:1) ,f , = r y + r z (cid:0) r y (2 r r y +2 r r y − rr y − rr y − r y − − r z ) − r r y +2 r r y − rr y − rr y − r y − r y − − r z r r y + − rr y − r y − r z − (cid:1) ,f , = r m r y r z r y + r z (cid:0) r (4 r − r y − r − − r z ) − rr y +2 r + r y − − r z ) + rr − (cid:1) ,f , = r m r y r z √ r y + r z (cid:0) r − − r z − rr y +2 r + r y r z − r y − − r z ) (cid:1) ,f , = r y + r z (cid:0) r y (2 r r y +6 r r y +4 r − rr y − r − r y +1)(1 − r z ) − r (1 − r z ) r − + r (1 − r z )(4 rr y +3 r − r − − r y (8 r r y +12 r − r r y − r r y − r +2 rr y +10 r + r y +1)( r − − r z ) + r r y − r r y − r r y − r +3 rr y + r + r y r − (cid:1) ,f , = r y + r z (cid:0) − r y (4 r r y +6 r − rr y − r − r y )1 − r z + r y (2 r r y +6 r r y +4 r − rr y − rr y − r +1)(1 − r z )
20 4 r r y − rr y − r y − r z (cid:1) . (24)Here, r = m b / ( m b + m cc