NLO QCD corrections to $J/\psi$ pair production in photon-photon collision

We calculate the next-to-leading order (NLO) quantum chromodynamics (QCD) correction to the exclusive process $\gamma+\gamma\to J/\psi+J/\psi$ 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 process will be hopefully to be seen by the Belle II detector within the next decade.


I. INTRODUCTION
The J/ψ meson has been proposed as an ideal laboratory to investigate both 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 the process dependent, but perturbative calculable short distance coefficients multiplied by the supposed 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 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 photon-photon scattering was studied extensively [2][3][4][4][5][6]. 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 the γ + γ → J/ψ + J/ψ process at the NLO accuracy.
The rest of the paper is organized as follows. In section II, we present the primary formulae employed in the calculation. In section III, we elucidate some technical details for the analytical calculation. In section IV, the corresponding numerical evaluation is performed. The last section is reserved for summary and conclusions.

II. FORMULATION
In this work, we investigate the photon-photon scattering in 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 Ee 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] f γ (x) = 1 where r = x xm(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 The total cross section can be expressed as the convolution between the cross section for γ + γ → J/ψ + J/ψ process and the photon distribution functions, where the dσ is calculated perturbatively up to the NLO level, The born level cross section and the NLO correction take the following forms: Hereŝ 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 andp are the momenta of heavy quark and antiquark respectively, P = p +p 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, 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 =p = P/2 and E = m c .
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.
In the calculation of NLO correction, the ultraviolet (UV) and infrared (IR) singularities are regularized by the dimensional regularization with D = 4 − 2ǫ. The IR singularities are canceled each other and the UV singularities are removed by renormalization procedure. The relevant renormalization constants include Z 2 , Z 3 , Z m and Z g , which corresponding to 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-massshell (OS) scheme, while others are defined in the modified minimal-subtraction (MS) scheme. The counter terms are 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 is 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 [21] is used to generate Feynman diagrams; FeynCalc [22] and FeynCalcFormLink [23] are used to handle the algebraic calculation; The package FIRE [24] 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 [25]. The overall phase space integrals are performed numerically by using the package CUBA [26].

IV. NUMERICAL RESULTS
In the numerical calculation, the general parameters are taken as α = 1 137.065 , m e = 0.511 MeV and m c = 1.5 GeV. The J/ψ radial wave function at the origin is extracted from J/ψ's leptonic width: By taking µ 0 = 2m c and Γ(J/ψ → e + e − ) = 5.55 keV [27], we obtain |R LO S (0)| 2 = 0.528 GeV 3 and |R NLO S (0)| 2 = 0.907 GeV 3 , which will be used in the LO and NLO calculation respectively. The two-loop formula for the running coupling constant is employed in the NLO calculation, in which, L = In the near future, it turns out the γ + γ → J/ψ + J/ψ process is hopeful to be observed by Belle II detector at the SuperKEKB electron-positron collider due to its high luminosity. The collider beam energies of the positron and electron are 4 GeV and 7 GeV respectively, which yields a center-of-mass energy of 10.6 GeV. The corresponding LO and NLO cross sections for γ + γ → J/ψ + J/ψ process with the WWA photon as the initial state are shown in Table I. The angular cut θ c of the WWA is set to be 83 mrad, the crossing angle between the positron and electron beams [28]. The transverse momentum cut p t ≥ 0.01 GeV is imposed on each individual J/ψ. In order to estimate the theoretical uncertainty, four reasonable choice for the renormalization scale [29], i.e. µ = 2m c , µ = 4m 2 c + p 2 t , µ = √ŝ /2 and µ = √ŝ are taken, where √ŝ is the center-of-mass energy of the photon-photon system. It can be seen that 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.
The instantaneous luminosity of the SuperKEKB collider may reach 8 ×  10.05) × 10 3 J/ψ pair events per year. In experiment, J/ψ can be reconstructed through its leptonic decays. According to Ref. [27], the branching ratio Br(J/ψ → l + l − (l = e, µ)) = 12%, then the number of candidates per year is 36 ∼ 145. Considering the signals are relatively clean, the theoretical estimation is hopefully testable in Belle II experiment.
The differential cross sections dσ/dcosθ at the SuperKEKB collider are shown in Fig.   2(a), where θ is the scattering angle. It can be seen that the events are more likely to be produced along the electron beam direction, which is one of the consequences of the asymmetry beams energy of the SuperKEKB collider. The differential cross sections dσ/dM J/ψJ/ψ are shown in Fig. 2(b), where M J/ψJ/ψ = √ŝ is invariant mass of the J/ψ pair. The invariant mass distribution arrives at the peak near threshold energy.
In the future e + e − colliders, like the Circular Electron Positron Collider (CEPC)  Table II. The angular cut θ c in WWA is set to be 32 mrad, as at the Large Electron-Positron Collider (LEP) [2].
Due to the LBS photon are more likely to be produced with large momentum fraction x, while the partonic cross sectionσ(γ + γ → J/ψ + J/ψ) arrives at the peak near  The kinematic cuts 2 ≤ p t ≤ 40 GeV and |y| < 2 are imposed on each single quarkonium. As further investigation, we also calculate the total cross sections for η c pair, Υ pair and η b pair production, as shown in Table II. It can be seen that the bottomonium pair production rates are 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. While for the η c pair production, the NLO cross sections are larger than that of the J/ψ pair production. However, due to the low reconstruction efficiency of η c , this process is also barely possible to be seen. In doing the numerical calculation for the bottomonium pair production, following parameters are used: m b

V. SUMMARY AND CONCLUSIONS
In this work, we investigate the J/ψ pair production in photon-photon fusion at the NLO QCD accuracy in the framework of NRQCD factorization formalism. The total cross section and the differential cross sections in bins of cosθ and M J/ψJ/ψ at the SuperKEBK collider are given. The total cross sections for J/ψ pair, η c pair, Υ pair and η b pair productions at the CEPC and the ILC are also estimated.
Numerical results shows that after including the radiative correction, the total cross sections are generally decreased, while their renormalization scale dependence is increased. This is a bit counterintuitive and needs a further investigation. The total cross section for J/ψ pair production via photon-photon fusion at the SuperKEBK collider is 0.100 ∼ 0.399 fb. And the corresponding event number is about (2.52 ∼ 10.05) × 10 3 per year. Considering J/ψ may be reconstructed via J/ψ → l + l − (l = e, µ) processes, the double J/ψ events may in the end reach 36 ∼ 145 a year. 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 NRQCD formalism to the NLO calculation of double quarkonium production poses another independent check on its validity at the next-to-leading order, especially with future experimental measurements.