Next-to-next-to-leading-order QCD corrections to double J/ψ production at the B factories

: In this paper, we study the next-to-next-to-leading-order (NNLO) QCD corrections for the process e + e − → J/ψ + J/ψ at the B factories. By including the NNLO corrections, the cross section turns negative due to the poor convergence of perturbative expansion. Consequently, to obtain a reasonable estimation for the cross section, the square of the amplitude up to NNLO is used. In addition, the contributions from the bottom quark and the light-by-light part, which are usually neglected, are also included. The ﬁnal cross section is obtained as 1 . 76 +2 . 42 − 1 . 66 fb at a center-of-mass energy of √ s = 10 . 58 GeV. Our result for total cross section and diﬀerential cross section could be compared with precise experimental measurement in future at the B factories.


Introduction
In quantum chromodynamics (QCD), the study of heavy quarkonium production assumes an important role in elucidating the interaction between quarks within two-body systems.In processes involving large momentum transfers, perturbative QCD is essential to estimate the theoretical results.In order to apply perturbative QCD to quarkonium production, various models have been introduced, including the color-evaporation model [1,2], the color-singlet model [3][4][5], and the nonrelativistic QCD (NRQCD) factorization formalism [6].The NRQCD factorization formalism, in particular, allows us to make consistent theoretical predictions and improve them order by order in the QCD coupling constant α s and the relative velocity of heavy quarks, denoted as v.
One of the most intriguing subjects within the domain of heavy quarkonium production and NRQCD is the phenomenon of double charmonium production in e + e − annihilation at the B factories.The experimental measurements for the processes e + e − → J/ψ + η c and e + e − → J/ψ + χ cJ at the B factories have been successfully performed by BELLE [7,8] and BABAR [9].However, in the case of the process e + e − → J/ψ + J/ψ, a clear signal is unable to be detected in BELLE's measurements, resulting in an upper limit of σ[e + e − → J/ψ + J/ψ] × B >2 < 9.1 fb [8].Here, B >2 signifies the branching ratio of final states involving more than 2 charged tracks.These experimental findings have sparked considerable theoretical endeavors, with the majority of investigations conducted within the framework of NRQCD factorization.Regarding the processes e + e − → J/ψ+η c and e + e − → J/ψ + χ cJ , their theoretical predictions have been computed up to two-loop level [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27].The results exhibit remarkable agreement with the experimental measurements, and further deepen our understanding of these interesting phenomena.
Regarding the process e + e − → J/ψ+J/ψ, the calculation at the NRQCD leading order (LO) provides a theoretical prediction for the total cross section about 8.7 fb [10].This value is even greater than the LO NRQCD prediction for the process e + e − → J/ψ + η c , and it has been updated to 6.65 fb shortly [28,29].With the aid of vector-dominance two-photon exchange model, the authors provide a total cross section about 2.38 fb by exclusively considering the photon fragmentation contribution [30].The non-fragmentation contribution has been investigated in Ref. [31] within the NRQCD factorization framework, and the authors find a significant destructive interference effect, resulting in a cross section prediction about 1.69 ± 0.35 fb.The next-to-leading-order (NLO) QCD correction for this process has been studied in Ref. [32], and it is found to be both negative and substantial.Including the NLO correction, the prediction shifts from the range of 7.4 ∼ 9.1 fb to −3.4 ∼ 2.3 fb.The results combining both NLO QCD and relativistic corrections are explored in Ref. [33].The authors find the fixed-order NRQCD prediction for the cross section is in the range of −12 ∼ −0.43 fb, which is negative and sensitive to the charm quark mass and renormalization scale.They also find that the predicted cross section will shift to the positive range of 1 ∼ 1.5 fb if the approach from Ref. [31] is used.
In this paper, the next-to-next-to-leading-order (NNLO) QCD corrections to e + e − → J/ψ+J/ψ are studied.It is found that the cross section becomes negative at NNLO, due to the poor convergence of perturbative expansion.In order to handle this, we use the square of the NNLO amplitude to obtain a reasonable result.In addition, the contributions from the bottom quark and light-by-light part, which are usually neglected, are also included.Recently, the NNLO corrections to this process are also studied in Ref. [34] where an improved NRQCD approach is proposed.
The remaining parts of the paper are organized as follows.In Section 2, we will provide relevant formulas and offer a concise overview of the calculation.Section 3 will be dedicated to presenting the numerical results and engaging in discussions.Section 4 will serve as the summary.

Cross section
Within the framework of NRQCD factorization, the cross section for e + (k 1 )e − (k 2 ) → J/ψ(p 1 ) + J/ψ(p 2 ) can be expressed as follows: where dσ represents the short-distance coefficients (SDCs), n 1,2 denotes all possible intermediate states, and O J/ψ (n 1,2 ) denotes the long-distance matrix elements (LDMEs).In the lowest-order nonrelativistic approximation, only the color-singlet state 3 S [1] 1 needs to be considered in the summation over n 1 and n 2 , hence we set As the LDME O J/ψ (n) incorporates nonperturbative hadronization effects, we initiate our calculation with the cross section of two on-shell (cc)-pairs with the quantum number 3 S [1] 1 .This cross section corresponds to the same SDCs as in e + e − → J/ψ + J/ψ, which can be expressed as: (2.2) Here, the symbol O (cc) [n] (n) is related to NRQCD bilinear operators and can be represented as: where the matrix element 0|χ † ǫ • σψ|cc(n) can be computed within the NRQCD framework [6,[35][36][37][38][39][40].On the other hand, the left-hand side of Eq.(2.2) can be directly calculated within perturbative QCD.Consequently, the SDCs dσ e + e − →(cc)[n]+(cc)[n] can be determined from Eq.(2.2).In combination with LDME O J/ψ (n) , we can obtain the cross section for e + e − → J/ψ + J/ψ as given in Eq.(2.1).
In the proceeding process, the e + e − pair initially annihilates into two virtual photons, which subsequently turn into two final states.The differential cross section can be conveniently expressed as follows: Here, the factor of 1/4 accounts for the spin average of the initial e + e − pair, 1/(2s) is the flux factor, s = (k 1 + k 2 ) 2 is the squared center-of-mass energy, A signifies the amplitudes for the process e + e − → (cc)[n] + (cc)[n], and dΦ 2 corresponds to the differential phase space for the two-body final state.

Calculation of the perturbative SDC
In this section, we will provide a concise overview of the calculation procedures.Firstly, we utilize the package FeynArts [41] to generate corresponding Feynman diagrams and amplitudes for e + e − → (cc)[n] + (cc)[n] at NNLO in α s .Secondly, we implement the package FeynCalc [42,43] to handle Lorentz index contraction and Dirac/SU (N c ) traces.Thirdly, we employ the package CalcLoop1 to decompose the Feynman amplitudes into 234 Feynman integral families, including the partial fraction decomposition of the linearly-dependent propagators in the integrals.Then the total NNLO amplitudes can be expressed by 87287 Feynman integrals which are further deduced into master integrals by using Kira [44] (a tool for integration-by-parts reduction).Finally, we utilize the package AMFlow [45][46][47][48][49] to calculate these master integrals.During the calculation, there are nearly 600 two-loop diagrams for e + e − → (cc)[n] + (cc) [n].Several representative Feynman diagrams up to two-loop order are illustrated in Figure 1.Here, p 3 = p 1 +q 1 2 , p 4 = p 1 −q 1 2 , p 5 = p 2 +q 2 2 , and p 6 = p 2 −q 2 2 , these momenta p 3,4,5,6 correspond to the charm and anti-charm quarks and satisfy the on-shell conditions p 2 3,4,5,6 = m 2 c .The momenta q 1 and q 2 represent the relative momenta between the quark and antiquark within the cc pairs.In the lowest-order nonrelativistic approximation, the relative momenta q 1 and q 2 are set to zero, which leads to p 2 1 = p 2 2 = (2m c ) 2 .We employ the conventional dimensional regularization approach with d = 4 − 2ǫ to regulate ultraviolet (UV) and infrared (IR) divergences.Feynman diagrams with a virtual gluon line connected to the quark pair in a meson exhibit Coulomb singularity, which manifest as power divergences in the IR limit of relative momentum.Such singularity can be addressed through cc wave function renormalization [17,50].In our calculation, we set the relative momenta q 1 and q 2 to zero before performing loop integration.The Coulomb divergence vanishes during the calculation with dimensional regularization.The UV divergence is resolved via renormalization.We employ on-shell (OS) renormalization for the heavy quark field and the heavy quark mass.The coupling constant α s is renormalized in the MS scheme.More explicitly, the amplitudes are renormalized according to: where the A il bare | i=0,1,2 represent the tree, one-loop, and two-loop bare amplitudes, respectively.The Z 2,c stands for the on-shell wave-function renormalization constant for the charm quark.The bare mass is renormalized as m Q,bare = Z m,Q m Q , with Z m,Q representing the on-shell mass renormalization constants for heavy quarks.The bare coupling constant is renormalized as: which corresponds to the MS scheme with n f active flavors.Here, µ R represents the renormalization scale, and Z MS αs stands for the renormalization constant of the coupling constant under the MS scheme.The renormalization constants Z MS αs up to two-loop level have been presented in Refs.[51][52][53][54][55]. Then the renormalized A(α s , m) can be obtained by expanding the right-hand side of Eq. (2.5) over renormalized quantities up to O(α 3 s ), i.e., Here, the A il | i=0,1,2 represent the tree, one-loop, and two-loop renormalized amplitudes, respectively.The loop integrals are computed with the measure µ 2ǫ R d d k/(2π) d , and the corresponding renormalization constants (Z 2,c , Z m,Q , and Z MS αs ) can be found in Refs.[56,57].Thus, the differential cross section can be obtained as: where κ = 1 − (16m 2 c )/s and θ is the angle between the J/ψ and the beam.To calculate the square of NNLO amplitude as shown in Eq. (2.8), we try to obtain the amplitude A nl | n=0,1,2 with the help of a complete basis space.These bases are constructed based on the Lorentz structures governing the process e + e − → (cc (2.9) Here, ρ 1 and ρ 2 represent the Lorentz indices of the two final states, u and v are the Dirac spinors of initial lepton pair in which spin notation is omitted, and the relationship p 1 = k 1 + k 2 − p 2 has been applied.It is worth noting that all the external particles in the process are colorless.Therefore, no color projection operator appears in the amplitude basis.This allows us to express the amplitude as follows: The coefficients c i are determined by where is the inner product of the amplitude A nl and the basis |e i , and G i,j = e j |e i signifies the inner product of the basis |e i and the basis |e j .It has been observed that the coefficients d 3 , d 6 , and d 8 are zero at the tree, one-loop, and twoloop levels, respectively.Then the products A ml A nl, * can be obtained using the following expression: where m and n take on values of 0, 1, and 2, respectively.However, there are remaining IR divergences in A 2l .This further makes A 2l A 0l, * , A 2l A 1l, * , and |A 2l | 2 divergent.On the other hand, as already known, O (cc)[n (n) also becomes IR divergent at two-loop level.In the leading order of v 2 and within the MS scheme, it can be expressed as which is derived from Refs.[6,[35][36][37][38][39][40].Here the term µ 2 Λ /µ 2 R −2ǫ arises from the evolution of the α s , from the factorization scale µ Λ to the renormalization scale µ R , since the correction is initially obtained at the scale µ Λ .The factor [e γ E /(4π)] −2ǫ is a consequence of the α s definition within the MS scheme, as given in Eq. (2.6).The divergence found in A 2l A nl, * is exactly same as this one, which renders the SDCs obtained from Eq.(2.2) without any divergences.Meanwhile this introduces an explicit logarithmic dependence on the NRQCD factorization scale µ Λ , which is on the order of ln(µ 2 Λ /m 2 ) in the SDCs.On the other hand, this µ Λ dependence can be completely cancelled by considering the µ Λ dependence of the LDMEs at fixed order.
The LDME O J/ψ (n) is often approximated as follows: where R J/ψ s (0) represents the wave function of J/ψ at the origin.Combining this approximation with Eqs.(2.1, 2.2, 2.8, 2.13), the differential cross section for e + e − → J/ψ + J/ψ can be expressed in the following form: Here, the terms f i | i=0,1,2,3,4 represent the SDCs at corresponding perturbative orders.It should be noted that the results for f 3 and f 4 are divergence free and gauge invariant, but incomplete, since only the contributions from A 2l A 1l, * and |A 2l | 2 are considered here.

Phenomenological results
For the numerical calculations, we use the following input parameters: Here, the bottom quark pole mass and running QCD coupling constant at the scale m Z are taken from Particle Data Group [58], and the QED coupling constant at the scale 2m c is taken from Refs.[15,34].m Z = 91.1876GeV is the Z boson mass and m c is the charm quark pole mass.We use the package RunDec3 [59] to evaluate the running QCD coupling constant α s (µ R ) at three-loop accuracy.
The value of R J/ψ s (0) can be extracted from the leptonic decay width at the two-loop level using the following formula [35,36,38,60,61]: where e c , C F , C A , T F , ζ 3 , and β0 = (11 − 2n f /3) /4 represent various constants and parameters.For the fermion-loop contributions to Γ J/ψ→e + e − , it is convenient to introduce Here n L = 3 counts the massless quarks, n H = 1 and n M = 1 label the contributions with closed massive quark loops with charm quark and bottom quark, respectively.The full expression of c b is shown in Ref. [61], whose numerical value is −0.394 in present case.By choosing Γ J/ψ→e + e − = 5.53KeV [58], µ Λ =1GeV as 1.810 GeV 3 at the two-loop level 2 .In the following, we consider the differential cross section with eleven sample points of | cos θ| ranging from 0 to 1.
In Tables 1, 2, and 3, we present the SDCs f i | i=0,1,2,3,4 defined in Eq. (2.15) for eleven different values of | cos θ|.It has been found that the numerical results of the SDCs f i | i=0,1,2,3,4 remain consistent whether we employ 10-digit or 20-digit precision for each Feynman integral family.In Table 1, there are five non-logarithmic terms in f 2 .The first three terms correspond to the vacuum polarisation and renormalisation contributions from the massless quark, charm quark and bottom quark, respectively.The fourth term corresponds to the light-by-light contributions from the massless quark (labeled as lbl).The fifth term represents all the other contributions, including the light-by-light contributions from charm quark and bottom quark.To check our calculation, we compare our numerical results at | cos θ| = 0.999 with the values presented in Table B.1 in the supplementary ma- 2 The leptonic decay width, denoted as Γ J/ψ→e + e − , is related to the decay constant f J/ψ through the relationship: As the same treatment as taken in Refs.[36,62,63], the leptonic decay widths at the n-loop level correspond to the calculation of the decay constant f J/ψ up to the n-loop level.Consequently, the predictions for the leptonic decay widths are essentially the square of the decay constant f J/ψ .If we derive the value of R J/ψ s (0) using the expression of the leptonic decay widths up to α 2 s -order, the numerical result for |R J/ψ s (0)| 2 = 5.528 GeV 3 is much bigger than the values estimated in various theoretical models [61,62,[64][65][66][67][68][69][70][71][72].A collection of the values of |R J/ψ s (0)| 2 are given in Ref. [63], which shows the |R (0)| 2 = 1.454GeV 3 in the Cornell potential model.This discrepancy is due to the poor convergence of perturbative expansion in the leptonic decay width Γ J/ψ→e + e − .However, this discrepancy can be mitigated by employing Eq. (3.2), which includes certain higher-order corrections in αs.In a similar vein, we study the square of NNLO amplitude for the process e + e − → J/ψ + J/ψ to obtain a reasonable theoretical prediction.
, and . In our results, we set n f = 5, n L = 3, n H = 1, n M = 1, and lbl = 1.terial of Ref. [34].From Table 1, we have πf 1 /f 0 = −10.78,which is exactly the same as that of Ref. [34].Taking the inputs n f = 4, n L = 3, n H = 1, n M = 0, and lbl = 0, we have π 2 f 2 /f 0 = −130.51− 22.46L µ − 51.18L µ Λ , which is consistent with that of Ref. [34].The light-by-light contributions are about 0.05% of the NNLO predictions, which are indeed small and negligible as claimed in Ref. [34].Furthermore, in Table 1, we observe that the numerical values of coefficients f 1 and f 2 are negative, resulting in negative corrections to LO predictions.In combination with α s ( √ s/2) = 0.1756, the theoretical prediction for the process e + e − → J/ψ + J/ψ at NNLO is unphysical due to poor convergence behavior.However, by incorporating two additional parts of higher-order corrections, as illustrated in Table 2 and Table 3, we are able to obtain a reasonable theoretical prediction.This is achieved by introducing positive corrections originating from f 3 and f 4 .
, and L µΛ = ln In Fig. 2, we illustrate the dependence of the renormalization scale µ R on the differential cross section for the process e + e − → J/ψ + J/ψ at various perturbative orders with | cos θ| = 0.800.It can be seen from the figure that the NNLO prediction exhibits the largest µ R dependence among all the results.It can also be found that the prediction based , and L µΛ = ln on the square of NNLO amplitude (denoted as S-NNLO3 ) displays a more pronounced sensitivity to variations in µ R when compared to both the NLO prediction and the prediction based on the square of NLO amplitude (denoted as S-NLO 3 ).The predictions based on S-NLO and S-NNLO closely approximate the LO prediction in the region µ R ∼ 6.5 GeV.This observation points to a notable convergence behavior in this particular region.In Fig. 3, we present the differential cross section for e + e − → J/ψ + J/ψ as a function

Summary
In summary, we have computed the NNLO QCD corrections for the production of J/ψ+J/ψ in e + e − annihilation at a center-of-mass energy of √ s = 10.58GeV.The numerical results of integrated cross section4 of e + e − → J/ψ +J/ψ with three typical renormalization scales µ R at different perturbative orders are shown in Table 4.It can be seen that the prediction at NNLO becomes negative due to the poor convergence of perturbative expansion.However, based on the square of NNLO amplitude, the theoretical prediction of the cross section becomes reasonable.The similar situation could be found in the case of NNLO result for Tables 1, 2, and 3.In other words, the trapezoidal rule is expressed as follows: x 2 x 1 f (x)dx = f (x 2 )+f (   and the second uncertainty is attributed to the method used to estimate the integrated cross section from the differential cross section.It should be pointed out that the center value of our prediction (1.76 fb) is lower than the value (2.13 fb) presented in Table .I of Ref. [34], and our result exhibits a more pronounced dependence on the choice of µ R .Furthermore, we find that the light-by-light contributions are indeed small and negligible as claimed in Ref. [34].Both of these theoretical predictions are lower than the upper limit of σ[e + e − → J/ψ + J/ψ] × B >2 < 9.1 fb.Our result for total cross section and differential cross section could be compared with precise experimental measurement in future at the B factories.

Figure 3 .
Figure 3.The differential cross section of e + e − → J/ψ + J/ψ as function of | cos θ| at various perturbative order, and the bands are obtained by varying the renormalization scale µ R within the range of [2m c , √ s].