NLO QCD corrections to exclusive quarkonium-pair production in photon–photon collision

We calculate the next-to-leading order (NLO) quantum chromodynamics (QCD) corrections to the exclusive processes γ+γ→Q+Q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma +\gamma \rightarrow {\mathcal {Q}}+{\mathcal {Q}}$$\end{document}, with Q=J/ψ,ηc,Υ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {Q}}=J/\psi ,\ \eta _c,\ \Upsilon $$\end{document}, or ηb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _b$$\end{document}, in the framework of non-relativistic QCD (NRQCD) factorization formalism. The cross sections at the SuperKEKB electron–positron collider, as well as at the future colliders, like the circular electron positron collider (CEPC) and the international linear collider (ILC), are evaluated. Numerical result indicates that the processes for J/ψ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J/\psi $$\end{document}-pair production and ηc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _c$$\end{document}-pair production are hopefully observable at the Belle II detector within the next decade.


Introduction
The J/ψ meson is thought of an ideal laboratory to investigate the perturbative and non-perturbative properties of quantum chromodynamics (QCD) since its discovery in 1974. The non-relativistic QCD (NRQCD) factorization formalism [1], which proposed by Bodwin, Braaten, and Lepage, provide a systematic framework for the theoretical study of quarkonium physics. In the NRQCD factorization formalism, the quarkonium production and decay rates can be factorized as process dependent, but perturbative calculable, short distance coefficients and universal long distance matrix elements (LDMEs). The relative importance between the LDMEs can be estimated by means of velocity scaling rules. In this way, the theoretical prediction takes the form of a double expansion in the strong coupling constant α s and the heavy quark relative velocity v.
Quarkonium production in photon-photon collision is an interesting topic to study, where the signals are relatively clean. At present, the inclusive J/ψ production via a e-mail: yanghao174@mails.ucas.ac.cn b e-mail: chenziqiang13@mails.ucas.ac.cn c e-mail: qiaocf@ucas.ac.cn (corresponding author) photon-photon scattering was studied extensively [2][3][4][5][6][7]. These researches indicate that the color-singlet (CS) sectors are inadequate to explain the experimental data [2,7], while the color-octet (CO) contributions based on different LDME sets diverse quite a lot, which even lead to opposite conclusions in some cases [2,6]. Besides the inclusive single J/ψ production, the exclusive J/ψ-pair production, where the CO contributions are suppressed by v 8 , also provides a test to the NRQCD factorization formalism at photon-photon colliders. In Ref. [8], the leading order (LO) calculation of exclusive J/ψ-pair production via photon-photon fusion was performed within the NRQCD factorization framework, while the J/ψ-pair diffractive production was investigated in the pomeron exchange scheme [9][10][11][12]. Considering the fact that the next-to-leading order (NLO) QCD corrections to physical processes in quarknoium energy regime are normally significant [13][14][15][16], in this work we calculate γ + γ → Q + Q processes, where Q = J/ψ, η c , ϒ, or η b , at the NLO accuracy.
The rest of the paper is organized as follows. In Sect. 2, we present the primary formulae employed in the calculation. In Sect. 3, we elucidate some technical details for the analytical calculation. In Sect. 4, the corresponding numerical evaluation is performed. The last section is reserved for summary and conclusions.

Formulation
In this work, we investigate the photon-photon scattering at e + e − colliders. The initial photons are assumed to be generated by the bremsstrahlung or the laser back scattering (LBS) effects. The energy spectrum of bremsstrahlung photon can be well formulated in Weizsacker-Williams approximation (WWA) [17] as where min , m e is the electron mass, x = E γ E e is the energy fraction of photon, √ s is the collision energy for e + e − collider, θ c is the maximum scattering angle of the electron or positron.
The energy spectrum of LBS photon is [18] where r = x x m (1−x) and N is the normalization factor: Here x m ≈ 4.83 [19] and the maximum energy fraction of LBS photon is restricted by 0 ≤ x ≤ x m 1+x m ≈ 0.83. The total cross sections can be expressed as the convolution of the γ + γ → Q + Q cross sections with the photon distribution functions, (4) where the dσ will be calculated perturbatively up to the NLO level, The Born level cross section and the NLO correction take the following forms: where 1 2 is a statistical factor for identical particles in the final state,ŝ is the center-of-mass energy square for the two photons, means sum (average) over the polarizations and colors of final (initial) state particles, dPS 2 denotes two-body phase space.
The standard form of quarkonium spin and color projection operator is adopted in our calculation [20]: Here, p Q and pQ are the momenta of heavy quark and antiquark respectively, P = p Q + pQ is the momentum of quarkonium, E = P 2 /4 is the energy of heavy (anti)quark in quarkonium rest frame, * S denotes the polarization vector of the spin-triplet state, 1 c represents the unit color matrix, and N c = 3 is the number of colors in QCD. At the leading order of the relative velocity expansion, it is legitimate to take p Q = pQ = P/2 and E = m Q . The projection operator for spin-singlet state can be readily obtained by replacing the * S with γ 5 in Eq. (7).

Analytical calculation
At the LO, there are twenty Feynman diagrams contributing to the γ + γ → Q + Q processes. Four of them are shown is Fig. 1. The LO calculation is straightforward and the result is pretty simple. We reproduce the Eq. (4) of Ref. [8], where the analytical expression for LO differential cross section was presented.
The typical Feynman diagrams for the NLO QCD corrections are shown in Fig. 2. According to their topologies, all diagrams can be classified into four groups: • Corrections to the tree-level Feynman diagrams, which are schematically shown in Fig. 2a-e. • Two-gluon exchange mechanism, as representatively shown in Fig. 2f-h. This mechanism is essential in vector Typical NLO Feynman diagrams for quarkonium-pair production via photon-photon fusion quarkonium-pair diffractive production, while forbidden in pseudo-scalar quarkonium-pair production due to the requirement of charge-parity (C-parity) conservation. • One of the quarkonia is produced through the two-gluon intermediate state, as schematically shown in Fig. 2ik. This topology is excluded in the calculation of vector quarkonium-pair production also due to the C-parity conservation. • Fermion box diagrams, as representatively shown in Fig. 2l.
Due to the fact that the C-parity selection rule excludes different topologies of vector and pseudo-scalar quarkonium-pair productions respectively, their NLO behaviors can be quite different.
In the calculation of NLO corrections, the ultraviolet (UV) and infrared (IR) singularities are regularized by the dimensional regularization with D = 4 − 2 . According to Ref. [1], the soft IR singularities are canceled with each other. Since the final charm quarks are massive, there is no collinear IR singularity in the hard-collinear region. The Coulomb IR singularities, which appears when a potential gluon is exchanged between the constituent quarks of a quarkonium, vanish in the dimensional regularization since we set p Q = pQ before calculating the Feynman integrals [21]. The UV singularities are removed by the standard renormalization procedure. The relevant renormalization constants include Z 2 , Z 3 , Z m and Z g , which correspond to the heavy quark field, gluon field, heavy quark mass and strong coupling constant, respectively. Among them, Z 2 and Z m are defined in the on-mass-shell (OS) scheme, while others are defined in the modified minimal-subtraction (MS) scheme. The counter terms are written as Here μ is the renormalization scale, γ E is the Euler's constant, β 0 = 11 3 C A − 4 3 T F n f is the one-loop coefficient of the QCD β-function, and n f denotes the active quark flavor numbers. The color factors C A = 3, C F = 4 3 and T F = 1 2 . In the calculation, the Mathematica package FeynArts [22] is used to generate Feynman diagrams; FeynCalc [23] and FeynCalcFormLink [24] are used to handle the algebraic calculation; The package FIRE [25] is employed to reduce all the one-loop integrals into typical master intergrals A 0 , B 0 , C 0 , D 0 , which can be numerically evaluated by the LoopTools [26]. The overall phase space integrals are performed numerically by using the package CUBA [27].

Numerical results
In the numerical calculation, the input parameters are taken as α = 1/137.065, m e = 0.511 MeV, m c = 1.5 GeV, Here, the J/ψ and ϒ radial wave functions at the origin are extracted from their leptonic widths with μ 0 = 2m Q , (J/ψ → e + e − ) = 5.55 keV, and (ϒ → e + e − ) = 1.34 keV [28]. Note, the LO and NLO extractions are employed in the LO and NLO calculation respectively. According to the heavy quark spin symmetry of NRQCD at the leading order in relative velocity expansion [1], here we take R η c (0) = R J/ψ (0) and R η b (0) = R ϒ (0).

The two-loop formula
for the running coupling constant is employed in the NLO calculation, in which, L = ln(μ 2 / 2 QCD ), β 0 = 11 , with n f = 4, QCD = 297 MeV for charmonium production, and n f = 5, QCD = 214 MeV for bottomonium production. For LO calculation, the one-loop formula for the running coupling constant is used.
According to our calculation, it turns out the charmoniumpair production process is hopeful to be observed by Belle II detector at the SuperKEKB electron-positron collider due to its high luminosity, where the beam energies of the positron and electron are 4 GeV and 7 GeV respectively, yielding a center-of-mass energy of 10.6 GeV. The corresponding LO and NLO cross sections for J/ψ-pair and η c -pair productions with the WWA photon as the initial state are shown in Table 1. The angular cut θ c of the WWA is set to be 83 mrad, the crossing angle between the positron and electron beams [29]. The transverse momentum cut p t ≥ 0.01 GeV is imposed on each individual quarkonium. In order to estimate the theoretical uncertainty, four reasonable choices for the renormalization scale [30], i.e. μ = 2m c , μ = 4m 2 c + p 2 t , μ = √ŝ /2 and μ = √ŝ are taken, where √ŝ is the centerof-mass energy of the photon-photon system. It can be seen that, for J/ψ-pair production, the total cross section are generally suppressed by the NLO corrections, and the theoretical uncertainty induced by renormalization scale is even larger at NLO than that at the LO, which seems a little counterintuitive and needs a further investigation. While for η c -pair production, after including the NLO corrections, the total cross sections are enhanced by a factor of about 2, and the theoretical uncertainty is reduced evidently.
The renormalization scale μ = 4m 2 c + p 2 t . The kinematic cuts 2 ≤ p t ≤ 40 GeV and |y| < 2 are imposed on each single quarkonium √ s (GeV) candidates per year would reach 36 ∼ 145 for J/ψ-pair production and 36-49 for η c -pair production. The differential cross sections dσ/dcosθ for J/ψ-pair and η c -pair productions at the SuperKEKB collider are shown in Fig. 3a, c respectively, where θ is the scattering angle. It can be seen that the events are more likely to be produced along the electron beam line, which is one of the consequences of the asymmetric beam energies of the SuperKEKB collider. The differential cross sections dσ/d M QQ for J/ψ-pair and η c -pair productions are shown in Fig. 3b, d respectively, where M QQ = √ŝ is invariant mass of the quarkonium-pair. The invariant mass distribution reaches the maximum around threshold energy.
Of the future high energy e + e − colliders, like the Circular electron positron collider (CEPC) [32] and the international linear collider (ILC) [33], the collision energy may reach 250 GeV or 500 GeV. And the LBS photon collision can be realized by imposing a laser beam to each e beam. Therefore, we investigate the J/ψ-pair and η c -pair productions under both WWA and LBS photon collisions with √ s = 250 GeV and √ s = 500 GeV. The corresponding LO and NLO total cross sections are presented in Table 2. The angular cut θ c in WWA is set to be 32 mrad, as at the Large Electron-positron Collider (LEP) [2]. Due to the fact, that the LBS photon are more likely to be produced with large momentum fraction x, while the partonic cross section σ (γ + γ → Q + Q) arrives at the peak around the threshold energy, the total cross sections for charmonium-pair production in the LBS photon case are only 1/10 to 1/2 the total cross sections in the WWA photon case. Given a typical collider luminosity 10 34 cm −2 s −1 ≈ 315 fb −1 year −1 , the event numbers of J/ψ-pair production and η c -pair production through WWA photon collision at √ s = 500 GeV are about 532 and 913 respectively. Assuming J/ψ is reconstructed through J/ψ → l + l − (l = e, μ), and η c is reconstructed through η c → KK π , several events might be observed per year in each case.
For completeness, we also calculate the total cross sections of ϒ-pair and η b -pair productions, as shown in Table 2. It can be seen that the bottomonium-pair production rates are two to three orders of magnitude smaller than the charmonium-pair production rates. Hence the bottomonium-pair production via photon-photon is invisible in the foreseeable future.

Summary and conclusions
In this work, we investigate the quarkonium-pair production in photon-photon fusion at the NLO QCD accuracy in the framework of NRQCD factorization formalism. The total cross sections and the differential cross sections in bins of cos θ and M QQ for charmonium-pair production at the SuperKEBK collider are given. The total cross sections for charmonium-pair and bottomonium-pair productions at the CEPC and the ILC are also estimated.
Numerical results shows that, for J/ψ-pair production, after including the radiative correction, the total cross sections are generally decreased and their renormalization scale dependence is increased. This is a bit counterintuitive and needs a further investigation. For η c -pair production, with the NLO corrections, the total cross sections are enhanced by a factor of about 2, and their renormalization scale dependence is reduced. The total cross sections for J/ψ-pair production and η c -pair production via photon-photon fusion at the SuperKEBK collider are 0.100-0.399 fb and 0.270-0.368 fb respectively. Considering J/ψ may be reconstructed via J/ψ → l + l − (l = e, μ) processes, η c may be reconstructed via η c → KK π , the double J/ψ and double η c events may in the end reach 36-145 and 36-49 a year, respectively. Since the signals are significant in experiment, the possibility of observing double quarkonium in future experiments still exists. Finally, it is remarkable that the successful application of the NRQCD factorization formalism to the NLO calculation of double quarkonium production poses another independent check on its validity at the next-to-leading order.