Analytical two-loop amplitudes of $e^{+} e^{-} \longrightarrow \boldsymbol{J} / \boldsymbol{\psi}+\boldsymbol{\eta}_c$ at $B$ factories

TL;DR

This paper presents an analytical NNLO calculation of two-loop amplitudes for $e^{+} e^{-} o J/\psi + \eta_c$ at B factories, improving theoretical predictions and aligning them with experimental data.

Contribution

The authors perform the first analytical NNLO calculation of these amplitudes, revealing the factorization structure and providing accurate cross section predictions.

Findings

Numerical predictions agree with experimental measurements.

Analysis of leading logarithms clarifies factorization in the process.

Extension to $Upsilon+eta_b$ production is discussed.

Abstract

In double charmonium production, a long-standing challenge is that the theoretical predictions are not consistent with the measurements at B factories. Within the NRQCD framework, the next-to-leading order (NLO) calculation has proved its power to cut down the discrepancy between theory and experiments. To further clarify this puzzle, we have performed the next-to-next-to-leading order (NNLO) calculation. The amplitude is obtained as an analytical asymptotic expansion in the ratio of the squared charm-quark mass over the squared center-of-mass energy, $m_c^2/s$. We investigate the origin of the leading logarithms by performing a region analysis, revealing the intricate factorization structure in this process. We provide numerical predictions on the total cross sections of $J/\psi+\eta_c$ production, which agree with the experimental results. Extension of our computation to $\Upsilon+\eta_b$ production is also discussed.

Tables

Table 1. Table 1 : Different input values taken for charm/bottom quarkonium production processes.

	$m_Q/\mathrm{GeV}$	$N_{lf}$	$N_{lp}$	$N_{lm}$	${\left|R_H\left(0\right)\right|}^2/{\mathrm{GeV}}^3$
$J/\psi +{\eta }_c$	1.5	3	1	2	0.978
${\rm Y}+{\eta }_b$	4.8	4	2	2	6.477

Table 2. Table 2 : The cross sections for e + ​ e − → J / ψ + η c e^{+}e^{-}\to J/\psi+\eta_{c} at different colliding energies. The central values are evaluated at μ R = s \mu_{R}=\sqrt{s} and m c = 1.5 m_{c}=1.5  GeV. The first uncertainty is estimated by varying μ R \mu_{R} from s / 2 \sqrt{s}/2 to 2 ​ s 2\sqrt{s} , and the second estimated by varying m c m_{c} from 1.3 1.3  GeV to 1.7 1.7  GeV.The factorization scale μ Λ = m c = 1.5 ​ GeV \mu_{\Lambda}=m_{c}=1.5~\text{GeV} .

$\sqrt{s}$ (GeV)	LO (fb)	NLO (fb)	NNLO (fb)
10.52	${3.93}_{-1.02-0.78}^{+1.70+0.73}$	${8.86}_{-1.98-2.11}^{+3.02+2.26}$	${12.29}_{-1.96-3.66}^{+2.07+4.50}$
10.58	${3.77}_{-0.98-0.74}^{+1.63+0.69}$	${8.53}_{-1.91-2.01}^{+2.90+2.15}$	${11.85}_{-1.90-3.49}^{+2.02+4.29}$
22.5	$\left({1.05}_{-0.23-0.03}^{+0.36+0.03}\right)\times {10}^{-2}$	$\left({2.87}_{-0.61-0.22}^{+0.92+0.24}\right)\times {10}^{-2}$	$\left({5.01}_{-0.96-0.62}^{+1.34+0.75}\right)\times {10}^{-2}$
91.2	$\left({9.83}_{-1.78-0.02}^{+2.46+0.02}\right)\times {10}^{-8}$	$\left({3.73}_{-0.73-0.16}^{+1.04+0.19}\right)\times {10}^{-7}$	$\left({8.48}_{-1.66-0.64}^{+2.34+0.75}\right)\times {10}^{-7}$
250	$\left({2.33}_{-0.37-0.00}^{+0.50+0.00}\right)\times {10}^{-11}$	$\left({10.87}_{-1.98-0.40}^{+2.73+0.47}\right)\times {10}^{-11}$	$\left({2.84}_{-0.53-0.18}^{+0.74+0.22}\right)\times {10}^{-10}$

Table 3. Table 3 : Cross sections for e + ​ e − → Υ + η b e^{+}e^{-}\to\Upsilon+\eta_{b} at different s \sqrt{s} . Central values are evaluated at μ R = s \mu_{R}=\sqrt{s} and m b = 4.8 m_{b}=4.8  GeV. The first uncertainty is from varying μ R \mu_{R} from s / 2 \sqrt{s}/2 to 2 ​ s 2\sqrt{s} , and the second from varying m b m_{b} from 4.6 4.6  GeV to 5.0 5.0  GeV.

$\sqrt{s}$ (GeV)	LO (fb)	NLO (fb)	NNLO (fb)
22.5	${0.073}_{-0.016-0.023}^{+0.025+0.025}$	${0.128}_{-0.020-0.043}^{+0.025+0.047}$	${0.115}_{+0.002-0.041}^{-0.023+0.048}$
91.2	$\left({4.06}_{-0.73-0.02}^{+1.01+0.02}\right)\times {10}^{-6}$	$\left({10.40}_{-1.72-0.20}^{+2.30+0.21}\right)\times {10}^{-6}$	$\left({16.28}_{-2.13-0.56}^{+2.51+0.59}\right)\times {10}^{-6}$
250	$\left({1.02}_{-0.16-0.00}^{+0.21+0.00}\right)\times {10}^{-9}$	$\left({3.31}_{-0.54-0.05}^{+0.72+0.05}\right)\times {10}^{-9}$	$\left({6.41}_{-0.96-0.16}^{+1.22+0.17}\right)\times {10}^{-9}$

Equations

$$\sigma (J/\psi +{\eta }_c(1S))\mathcal{B}(2)=56\pm 8\pm 4\mathrm{fb},$$

$$\sigma (J/\psi +{\eta }_c(2S))\mathcal{B}(2)=65\pm 0\pm 4\mathrm{fb},$$

$$\sigma \left(J/\psi +{\eta }_c(1S)\right)\mathcal{B}(2)=76\pm 8_{-2.1}^{+1.5}\mathrm{fb},$$

$$\sigma \left(J/\psi +{\eta }_c(2S)\right)\mathcal{B}(2)=64\pm 7_{-3.0}^{+2.4}\mathrm{fb},$$

$$e^{+}(k_1)+e^{-}(k_2)\to J/\psi (P_1)+{\eta }_c(P_2).$$

$$k_1^2=k_2^2=0,P_1^2=P_2^2=4m_c^2,(k_1+k_2{)}^2=s,$$

$${\mathcal{A}}^{\mu }(s)=⟨J/\psi (P_1,\lambda )+{\eta }_c(P_2)\left|J^{\mu }\right|0⟩$$

$$=$$

$$\times \sqrt{M_{J/\psi }{M}_{{\eta }_c}}{R}_{J/\psi }(0)R_{{\eta }_c}(0),$$

$$\sigma =\frac{\alpha }{12s^2}\sqrt{\frac{s-16m_c^2}{s}}{\mathcal{A}}^{\mu }{\mathcal{A}}_{\mu }^{\ast }.$$

$$h^{(0)}(r)=-i\frac{512\left({\mathrm{N}}_c^2-1\right)}{3}{r}^2.$$

$$\mathcal{I}(r)=\sum _{k,n}{c}_{k,n}(ϵ){\ln }^k(r)r^n.$$

$${\tilde{\mathcal{A}}}^{\mu }({\alpha }_s,m_Q)=Z_{2,Q}^{-2}\left({\mathcal{A}}_0^{\mu (0)}+{\mathcal{A}}_0^{\mu (1)}+{\mathcal{A}}_0^{\mu (2)}\right),$$

$$m_{Q,0}=Z_{m,Q}{m}_Q,{\alpha }_{s,0}=(4\pi e^{-{\gamma }_E}{)}^{-ϵ}{\mu }_R^{2ϵ}{Z}_{{\alpha }_s}^{\overline{\mathrm{MS}}}{\alpha }_s({\mu }_R),$$

$${\tilde{\mathcal{A}}}^{\mu }=-\frac{1}{2ϵ}\frac{{\alpha }_s^2}{{\pi }^2}({\gamma }_{J/\psi }+{\gamma }_{{\eta }_c}){\mathcal{A}}^{\mu (0)}+\mathcal{O}({ϵ}^0),$$

$${\gamma }_{J/\psi }$$

$${\gamma }_{{\eta }_c}$$

$${\bar{\mathcal{A}}}_{\text{finite}}^{\mu (2)}={\tilde{\mathcal{A}}}^{\mu (2)}+\frac{1}{2ϵ}\frac{{\alpha }_s^2}{{\pi }^2}({\gamma }_{J/\psi }+{\gamma }_{{\eta }_c}){\mathcal{A}}^{\mu (0)}{\left(\frac{{\mu }_{\Lambda }^2}{{\mu }_R^2}\right)}^{-2ϵ}$$

$$h^{(0)}(r)$$

$$h^{(1)}(r)$$

$$+\mathcal{O}(r^3),$$

$$h^{(2)}(r)$$

$$+L_{(2,2)}^{(2)}{\ln }^2(r)+L_{(2,1)}^{(2)}\ln (r)+L_{(2,0)}^{(2)}]$$

$$+r^{5/2}{L}_{(5/2,0)}^{(2)}+\mathcal{O}(r^3),$$

$$L_{(2,2)}^{(1)}=$$

$$L_{(2,1)}^{(1)}=$$

$$-\frac{4i({\mathrm{N}}_c^2-1)}{9{\pi }^2{\mathrm{N}}_c}[(31+138\ln 2){\mathrm{N}}_c^2$$

$$-28{\mathrm{N}}_c-258\ln 2+57],$$

$$L_{(2,4)}^{(2)}=-\frac{i({\mathrm{N}}_c^2-1{)}^2(7{\mathrm{N}}_c^2-19)}{144{\pi }^4{\mathrm{N}}_c^2},$$

$$L_{(2,3)}^{(2)}=$$

$$+\frac{({\mathrm{N}}_c^2-1{)}^2\left(7{\mathrm{N}}_c^2-19\right)}{36{\pi }^3{\mathrm{N}}_c^2}$$

$$-\frac{i\left({\mathrm{N}}_c^2-1\right)}{54{\pi }^4{\mathrm{N}}_c^2}[{\mathrm{N}}_c^4(79+87\ln 2)-31{\mathrm{N}}_c^3$$

$$-2{\mathrm{N}}_c^2(40+141\ln 2)+{\mathrm{N}}_c(59-24\ln 2)$$

$$-99+189\ln 2].$$

$$L_{(5/2,0)}^{(2)}=-\frac{8i({\mathrm{N}}_c^2-1)(7{\mathrm{N}}_c^2-2)}{9{\pi }^2{\mathrm{N}}_c}.$$

$$h^{(0,2)}(r)$$

$$h^{(1,2)}(r)$$

$$h^{(2,2)}(r)$$

$$+L_{(2,2)}^{(2)}{\ln }^2(r)+L_{(2,1)}^{(2)}\ln (r)+L_{(2,0)}^{(2)}.$$

$$\mathrm{Im}\{h^{(0,2)}\}=$$

Full text

I Introduction

Back in the beginning of this century, exclusive double charmonium production in electron-positron annihilation at B factories has drawn great theoretical interest Abe et al. (2002, 2004); Aubert et al. (2005). Since these sorts of hard exclusive reactions serve as a novel and fertile ground for exploring the interplay between perturbative and non-perturbative aspects in the heavy flavor physics of QCD, many studies of double quarkonium production including $e^{+}e^{-}\rightarrow J/\psi+H$ have been performed by several research groups Braaten and Lee (2003a); Ma and Si (2004); Bondar and Chernyak (2005); Isgur and Smith (1984), where $H$ denotes a C-even charmonium that recoils against $J/\psi$. Recently, the Belle II data-taking is in progress Yin et al. (2023), which allows to probe the heavy flavor physics with unprecedented precision Altmannshofer et al. (2019). Future high-luminosity lepton colliders such as the ILC, CEPC, and FCC-ee are expected to provide complementary opportunities to study double charmonium production Dong et al. (2018); Durán et al. (2025); Abada et al. (2019). Thus, a reliable theoretical prediction is highly demanded.

Pioneering theoretical effort, including leading order (LO) NRQCD, light-cone or general pQCD approaches Braaten and Lee (2003a); Ma and Si (2004); Bondar and Chernyak (2005); Isgur and Smith (1984), exhibited unsatisfactory large uncertainties and discrepancies with experimental data. Especially, the LO NRQCD predictions Braaten and Lee (2003a) of $e^{+}e^{-}\rightarrow J/\psi+\eta_{c}$ is much smaller than the measured cross section Abe et al. (2002); Aubert et al. (2005). Here we list the experimental results:

[TABLE]

from Belle, and

[TABLE]

from BaBar, where $\mathcal{B}(n)$ denotes the branching fraction for $\eta_{c}$ decaying into more than $n$ charged tracks. In principle, the NRQCD theory Bodwin et al. (1995) provides an effective framework that allows us to systematically deal with the higher-order corrections, which may be helpful to understand this discrepancy.

The next-to-leading order (NLO) perturbative calculations within the NRQCD framework were computed in Ref. Zhang et al. (2006), which, for the first time, reaches the lower uncertainty bound of the experimental measurements. This computation demonstrates the importance of high-order corrections in the $e^{+}e^{-}\rightarrow J/\psi+\eta_{c}$ process. The NLO calculation was later verified independently by Gong and Wang (2008). Subsequent developments addressed the relativistic corrections for the process up to $\mathcal{O}(v^{2})$ Braaten and Lee (2003b); He et al. (2007), where $v$ represents the characteristic velocity of the constituent charm quarks in the quarkonium rest frame. While the non-perturbative NRQCD matrix elements inherently carry significant theoretical uncertainties, comprehensive studies demonstrated that combining both NLO perturbative corrections and partial resummation of relativistic effects Bodwin et al. (2008) could yield predictions that were broadly consistent with B-factory measurements, albeit with substantial residual uncertainties. Further refinement came with the investigation of joint perturbative and relativistic corrections at $\mathcal{O}(\alpha_{s}v^{2})$, which were found to provide modest but meaningful enhancements to the NRQCD theoretical predictions Dong et al. (2012).

The next-to-next-to-leading order (NNLO) radiative corrections were numerically computed in Ref. Feng et al. (2019) and further confirmed in Ref. Huang et al. (2023), which were found to be substantial. Although the inclusion of both radiative and relativistic effects brings the theoretical prediction into consistent with experimental measurements, significant uncertainties persist, primarily due to the choice of renormalization scale. It is therefore essential to investigate the origin of the poor perturbative convergence in this process. At NLO, the dominant contribution comes from the leading logarithm, as analyzed in Refs. Jia et al. (2011); Bodwin et al. (2014). However, the significance of logarithmic contributions at NNLO remains unclear.

In this paper, we provide a detailed computation of the NNLO corrections of the process $e^{+}e^{-}\rightarrow J/\psi+\eta_{c}$, utilizing state-of-the-art techniques. Our results are expressed as an asymptotic expansion in the variable $r=m_{c}^{2}/s$, enabling us to extract all logarithmic terms of $r$. The concise leading logarithmic terms are shown explicitly in the main text, while more extensive subleading terms are provided in an attached supplementary file. Furthermore, we explore the origin and structure of the leading logarithmic contributions and summarize the resulting theoretical insights. Using our analytical results, we provide physical predictions for the $J/\psi+\eta_{c}$ production cross section at various center-of-mass energies. As a nontrivial byproduct, our formalism also allows for a convenient prediction of the $\Upsilon+\eta_{b}$ cross section, which may be tested in future experiments. We additionally investigate the impact of input parameters and scale choices on the phenomenological outcomes and assess the convergence behavior of the perturbative expansion. The paper concludes with a summary of our main results.

II Framework and Methods

We consider the two-loop QCD correction to the following process

[TABLE]

At the lowest order in velocity expansion within NRQCD, the scalar products of external momenta satisfy

[TABLE]

where $\sqrt{s}$ is the center-of-mass energy and $m_{c}$ denotes the charm-quark mass. It has been proved that the cross section of $e^{+}e^{-}\rightarrow J/\psi+\eta_{c}$ can be factorized into the leptonic and hadronic parts in Ref. Sun (2021).

It is convenient to define the amplitudes of the process $\gamma^{\ast}\rightarrow J/\psi+\eta_{c}$ as

[TABLE]

with $J^{\mu}$ being the electromagnetic current and $h^{(i)}(r)$ being dimensionless expansion terms in $\alpha_{s}$. The quarkonium mass can be taken as $M_{H}^{2}\approx 4m_{c}^{2}$ and $R_{H}(0)$ denotes the wave function at the origin. Typical tree-level and one-loop Feynman diagrams are shown in Fig. 1.

With the amplitudes at hand, the total cross section for $e^{+}e^{-}\rightarrow J/\psi+\eta_{c}$ can be written as

[TABLE]

At the LO, we have

[TABLE]

In the following, we will compute $h^{(1)}(r)$ and $h^{(2)}(r)$.

We generate the Feynman diagrams for the process with FeynArts Nogueira (1993) and QGRAF Hahn (2001), with some representative two-loop diagrams shown in Fig. 2. Our calculation employs the dimensional regularization with $D=4-2\epsilon$. The amplitudes are computed using the CalcLoop package cal , which performs the Dirac algebra, Lorentz contractions, as well as $SU(N_{c})$ color algebra operations. Through tensor decomposition, we express the amplitudes as linear combinations of scalar Feynman integrals, which are classified into 167 distinct integral families. For each family, we utilize the packages Blade Guan et al. (2025), FiniteFlow Peraro (2016, 2019) and LiteRed Lee (2012) to perform the integration-by-part (IBP) Tkachov (1981); Chetyrkin and Tkachov (1981); Laporta (2000) reduction to express these integrals as linear combinations of a set of so-called master integrals. In our computation, we assign $\epsilon$ small numerical values and fit the $\epsilon$ dependence at the final stage as proposed in Refs. Liu and Ma (2022a, b), which avoids the manipulation of $\epsilon$ in the intermediate step and reduces the computational time.

To evaluate the master integrals, we use the method of numeric differential equations Kotikov (1991). We construct the differential equations with respect to the dimensionless variable $r=m_{c}^{2}/s$. The boundary conditions at $r=1/100$ are obtained using the auxiliary mass flow method Liu et al. (2018); Liu and Ma (2022c, a), implemented in AMFlow Liu and Ma (2022b). We note that spurious or non-physical singularities can cause numerical instability, which can be resolved by increasing the working precision in our numerical computation. Feynman integrals exhibit the following power series expansion structure Caffo et al. (2009); Czakon (2008); Lee (2015); Moriello (2020); Hidding (2021); Armadillo et al. (2023):

[TABLE]

We utilize the package DESolver implemented in AMFlow to solve the differential equation as a power series.

The precision of the resulting asymptotic expansion coefficients $c_{k,n}$ depends on both the accuracy of the boundary conditions and the truncation order of $r$. We truncate the expansion to at most $O(r^{2000})$ to ensure that the leading expansion coefficients have at least 600 digits of precision. The high precision results enable us to reconstruct the analytic form of the leading coefficients using the PSLQ algorithm Ferguson and Bailey (1992). The basis we use in this work for the PSLQ reconstruction of finite terms is the combination of three basic elements $\mathrm{Li}_{m}(-2),\mathrm{Li}_{m}(-1),\mathrm{Li}_{m}(2)$ for $m\leq 4$. With the analytic form of the leading coefficients as boundary conditions, we can then solve the differential equations again to obtain the analytic form for all other coefficients. Finally, we express the amplitudes as a series expansion in $r$ with fully analytical coefficients, enabling exact symbolic manipulation and further theoretical analysis.

III Renormalization and factorization

The amplitude, before factorization matching, can be written as

[TABLE]

where $Z_{2,Q}$ is the field renormalization constant for the heavy quarks in the final state, and $\mathcal{A}_{0}^{\mu(i)}$ are amputated bare amplitudes which in general depend on bare $\alpha_{s,0}$ and $m_{Q,0}$. The bare and renormalized quantities are related as

[TABLE]

where $\mu_{R}$ is the renormalization scale. The renormalization constants can be found in Refs. Czakon et al. (2008); Czakon (2008); Bärnreuther et al. (2014) and references therein.

After performing the UV renormalization, we are left with a single $\epsilon$ pole in the amplitude, which satisfies exactly the following factorization relation

[TABLE]

where the anomalous dimensions are defined as Czarnecki and Melnikov (2001); Kniehl et al. (2006); Beneke et al. (1998); Czarnecki and Melnikov (1998)

[TABLE]

The exact subtraction of the poles from the amplitudes provides a robust examination of our calculation. The finite amplitude after factorization matching reads:

[TABLE]

with $\mu_{\Lambda}$ being the factorization scale. Now, we have obtained the renormalized two-loop amplitude $\bar{\mathcal{A}}_{\text{finite}}^{\mu(2)}$, which is free of divergences. Unless otherwise specified, we simply use $\mathcal{A}^{\mu(i)}$ to indicate the finite amplitude in the following text. The corresponding $h^{(i)}(r)$ functions likewise have no divergences.

IV Analytical Results

We express $h^{(i)}(r)$ explicitly as following:

[TABLE]

where higher-power terms in $r$ expansion are ignored since they are irrelevant for actual phenomenological study.

The LO result is given in Eq. (7). The logarithmic terms at NLO read

[TABLE]

[TABLE]

where $N_{lf}$ denotes the number of light fermions. These results are found to be consistent with Refs. Jia et al. (2011); Gong and Wang (2008); Bodwin et al. (2014). The leading and next-to-leading logarithmic terms at NNLO read

[TABLE]

[TABLE]

Besides, the asymptotic expansion at two-loop order has an unusual non-integer power term $r^{5/2}$, with its coefficient given by

[TABLE]

Other coefficients are too lengthy to be given here, and they are provided in the supplementary materials.

We define $h^{(i,2)}(r)$ as the coefficients at leading power, given by

[TABLE]

By setting $\mu_{R}=\sqrt{s}$, $\mu_{\Lambda}=m_{c}$ and $N_{lf}=3$, we have the numerical results

[TABLE]

where only the imaginary part is given. With the above results, we can analyze the relative contribution from each logarithmic term at different energy points. As shown in Fig. 3, the $\ln^{2}(r)$ term dominates the NLO results. At NNLO, the $\ln^{3}(r)$ term dominates across most of the energy region of phenomenological interest, while the $\ln^{4}(r)$ term also contributes substantially at relatively higher energies. The observation implies that resummation of both the leading and next-to-leading logarithmic terms to all orders of $\alpha_{s}$ should be necessary.

V The structure of large logarithms

The large logarithms shown in the last section not only give a major contribution to the cross section, but also possess intriguing features that deserve detailed exploration. The leading logarithms at NLO have been analysed in Ref. Bodwin et al. (2014). It was found that the double logarithms originate solely from the soft fermion contribution while those from soft gluon contribution cancel out. This kind of logarithm has received attention recently and has been discussed in $H\to\gamma\gamma$ Liu and Neubert (2020); Wang (2019); Hou et al. (2025), $Hgg$ form factor Liu and Penin (2018); Liu et al. (2023), gluon thrust Beneke et al. (2022), and $e\mu$ backward scattering Bell et al. (2022), $B_{c}\to\eta_{c}$ form factors Bell et al. (2024), etc. Below we will provide an analytic calculation of the leading logarithms at both NLO and NNLO by using a method developed in Refs. Wang (2019); Hou et al. (2025). This would not only serve as a nontrivial check of the above result from differential equations, but also may help to develop a formalism for the resummation of such logarithms to all orders.

The large logarithms appear for the first time at one-loop. The double logarithms can arise when the gluon momentum becomes soft and (anti-)collinear; see Ref. ter Hoeve et al. (2024) for the analyses of the massive form factor with the method of regions Smirnov (1997); Beneke and Smirnov (1998). This contribution is proportional to the color charge of the external collinear particle and thus vanishes after summing over all diagrams, since the mesons are color neutral. The other kind of region that contributes to the double logarithms involves the soft quark. The relevant Feynman diagrams are shown in Fig. 4.

Take diagram 4(d) as an example. The loop momentum of the quark is denoted by $l$. The collinear and anti-collinear directions are labeled by $n^{\mu}=(1,0,0,1),$ and $\bar{n}^{\mu}=(1,0,0,-1)$, respectively. Any momentum $k$ can be represented by $(k^{+},k^{-},k_{\perp}^{\mu})$ with $k^{+}=n\cdot k$ and $k^{-}=\bar{n}\cdot k$. The external charm (anti-)quarks have collinear $p\sim(r,1,0)\sqrt{s}$ or anti-collinear $\bar{p}\sim(1,r,0)\sqrt{s}$ momenta. The double logarithms receive contributions from four regions: the hard region $l\sim(1,1,1)\sqrt{s}$, the collinear region $l\sim(r,1,\sqrt{r})\sqrt{s}$, the anti-collinear region $l\sim(1,r,\sqrt{r})\sqrt{s}$ and the soft region $l\sim(\sqrt{r},\sqrt{r},\sqrt{r})\sqrt{s}$. The soft region gives the amplitude (with the LO amplitude split)

[TABLE]

where we have neglected terms that do not contribute to the double logarithm. The Feynman prescription for the three denominators is given by $+i\varepsilon,+i\varepsilon,-i\varepsilon$, respectively. Naive integration leads to endpoint divergences even with a dimensional regulator. Following the method in Ref. Hou et al. (2025), we add auxiliary $\delta$-regulators in the full denominators, i.e., $(l+p)^{2}\to(l+p)^{2}+\delta_{1}s$ and $(l-\bar{p})^{2}\to(l-\bar{p})^{2}+\delta_{2}s$. Consequently, the above denominators in the soft region are modified as $l^{+}\to l^{+}+2\delta_{1}\sqrt{s},l^{-}\to l^{-}-2\delta_{2}\sqrt{s}$. Integrating over $l^{-}$ using the residue theorem, we then obtain

[TABLE]

The endpoint divergences at $l^{+}=0$ and $l^{+}=\infty$ are now regulated by $\delta_{1}$ and $\delta_{2}$, respectively. The integration of $l_{T}^{2}$ generates an ultraviolet divergence in $1/\epsilon$. The collinear and anti-collinear regions have the same structure. All the ultraviolet divergences would cancel with the infrared divergences in the hard region. Given this insight, we split the $l_{T}^{2}$ integration into two parts, $(0,s)$ and $(s,\infty)$. In the second part, $l_{T}$ is considered as a hard scale, and thus the small scale $m_{c}$ can be dropped. As a result, no logarithms would appear in this part. The first part contains no dimensional divergences and is given by

[TABLE]

where we have kept only the double logarithms. Similarly, we obtain the results in the collinear and anti-collinear regions,

[TABLE]

Note that the expansion of $\mathcal{A}_{c}$ in the soft (or anti-collinear) region reproduces $\mathcal{A}_{s}$. Therefore, a subtraction needs to be performed to avoid double-counting Manohar and Stewart (2007). Taking the sum of the above contributions yields

[TABLE]

which agrees with the expansion of the full one-loop amplitude. It can be seen that all the $\delta$-regulators cancel out. Another important feature of this method is that the collinear and anti-collinear results are proportional to $\ln(\delta_{i})$ and thus can be taken to be vanishing if setting $\delta_{i}=1$. This indicates a shortcut to obtain the leading logarithms. One needs to calculate only the soft region contribution analytically.

The other diagrams in Fig. 4 can be computed similarly. Based on these analyses of regions, we show the reduced diagrams contributing to double logarithms in Fig. 5. The whole perturbative amplitude can be considered as a convolution of the hard function, the radiative jet function, and the soft function, which is interfaced with projectors. The hard function represents the contribution of all hard propagators with $\mathcal{O}(1)$ off-shellness. It can be considered as the Wilson coefficients of the effective operators, consisting of the building blocks $W_{hc}^{\dagger}\xi_{hc}$ and $W_{hc}^{\dagger}[iD_{\perp}^{\mu}W_{hc}]$ Beneke et al. (2018a, b). The radiative jet function $J_{c,q}$ describes the transition of a hard collinear gluon to a collinear anti-charm-quark and a soft charm-quark, denoted by its subscripts $c$ and $q$, respectively. The other radiative jet function $J_{c\bar{c},q}$ describes the transition of a hard collinear charm-quark to a collinear (anti-)charm-quark pair and a soft charm-quark. The soft function $S$ incorporates the effects of all soft particles, including the soft Wilson lines which are not shown explicitly in the figures. The external collinear (anti-)charm-quark pair forms a state of spin-singlet (triplet) via the projectors Bodwin et al. (2014)

[TABLE]

where we have performed the power expansion in $r$ and the first two terms are given by

[TABLE]

In Fig. 5(a), the radiative jet and soft functions are of order $r^{-1/2}$ and $r^{1/2}$, respectively, and thus the soft loop momentum convolution gives a result of order $r^{1/2}$. The two leading power projectors are of order $r^{-1}$. So the amplitude of this diagram is of order $r^{-1/2}$, same as the LO amplitude. It is also interesting to note that the Dirac structures of $\Pi_{1}^{(0)}$ and $\Pi_{3}^{(0)}$ require a helicity flip along the fermion line, which can only occur in the soft function since the radiative jet functions are at leading power. In Fig. 5(b), the projector $\Pi_{3}^{(1)}$ is power suppressed. However, $J_{c\bar{c},q}^{(0)}$ is power enhanced by $r^{-1/2}$ compared to $J_{c,q}^{(0)}$. Thus, the amplitude is at the same power as Fig. 5(a). Fig. 5(c) contains two types of contributions. The first type, with $i=0$ and $j=1$, is analogous to Fig. 5(b) with the collinear and anti-collinear directions exchanged. The other type, with $i=1$ and $j=0$, is similar to Fig. 5(c), with the replacement of $\bar{J}_{c,q}^{(0)}$ by $\bar{J}_{c\bar{c},q}^{(1)}$, which is of the same power. The power counting of Fig. 5(d) can be understood following the same logic.

At the two-loop level, there are additional types of contributions. Sample reduced diagrams are shown in Fig. 6. Fig. 6(a) generalizes Fig. 5(d) to include the case with leading power projector $\Pi_{3}^{(0)}$ but with subleading radiative jet function $J_{c\bar{c},q}^{(1)}$. Fig. 6(b) introduces a new type of radiative jet function, $J_{c\bar{c},qg}$, that emits both a soft quark and a soft physical gluon. It can be considered as obtained by inserting in $J_{c\bar{c},q}$ a subleading power interaction ($\mathcal{O}(r^{1/4})$) involving a soft gluon. The suppression induced by $J_{c\bar{c},qg}$ is balanced by the enhancement of the leading power projector $\Pi_{3}^{(0)}$, compared to Fig. 5(d). The new factorization topologies appearing at higher loops represent a challenge to the development of a resummation formalism in this process.

We calculated most of the two-loop diagrams contributing to the leading logarithms with the method described above. Although the soft quark momentum mode is fixed, the additional gluon loop momentum can have hard, (anti-)hard-collinear, (anti-)collinear, (anti-)soft-collinear (scaling as $(r^{3/2},r^{1/2},r)$), and soft modes. The (anti-)soft-collinear mode does not appear in the reduced diagrams because its contribution cancels out in the amplitude due to the same reason as the soft gluon at one loop. We performed the calculation loop-by-loop. The additional loop was calculated analytically using the Feynman parameterization. Thanks to the $\delta$-regulators and the cut on $l_{T}^{2}$ applied to the remaining soft quark loop, the result can be safely expanded in $\epsilon$. Then the integration over $l$ was carried out using the residue theorem. We investigated 80 Feynman diagrams, finding agreement with the direct computation. It is crucial to compare the subset of diagrams in which the leading logarithms induced by pure soft gluons cancel.

However, there exist a few diagrams that cannot be computed using the $\delta$-regulators. These diagrams contain the tent-type subdiagram, which has also appeared in the muon-electron backward scattering Bell et al. (2022). The structure of the involved two-dimensional endpoint divergences requires advanced applications of the $\delta$-regulators (or other regulators). This interesting problem is left to future work. Nevertheless, the region analyses expose the intricate factorization structure hidden behind the leading logarithms and may shed light on the resummation of these logarithms to all orders.

VI Phenomenology

To facilitate comparison between theoretical predictions and experimental results, we summarize the input parameters used in this work in Tab. 1. The masses of charm quark and bottom quark are taken as $m_{c}=1.5$ GeV and $m_{b}=4.8$ GeV, respectively, while other light quarks are considered as massless in this paper. $N_{lf}$ denotes the number of active flavors of quarks, while $N_{lp}$ and $N_{lm}$ denote the number of active positive-charged and negative-charged quarks, respectively. $|R_{H}(0)|^{2}$ is the long-distance matrix element for the quarkonium $H$ in the NRQCD framework, which was taken from Refs. Zhang et al. (2006); Eichten and Quigg (1995); Brambilla et al. (2023). At the lowest order in velocity expansion, we have $|R_{J/\psi}(0)|^{2}=|R_{\eta_{c}}(0)|^{2}$ and $|R_{\Upsilon}(0)|^{2}=|R_{\eta_{b}}(0)|^{2}$.

The QCD coupling constant $\alpha_{s}(\mu_{R})$ is evaluated using our in-house program with $m_{c},m_{b}$ and $\mu_{R}$ as input, yielding, for instance, $\alpha_{s}(91.2~\text{GeV})=0.118$ and $\alpha_{s}(10.58~\text{GeV})=0.177$. Throughout this paper, the factorization scale $\mu_{\Lambda}$ is set to be $m_{Q}$. The QED coupling constant is fixed at $\alpha=1/137$. We also vary the quark mass by $\pm 0.2~\text{GeV}$ to analyze its uncertainty.

VI.1 $J/\psi+\eta_{c}$ production

We present in Tab. 2 the total cross sections at various collider energies. The cross section at $\sqrt{s}=10.58$ GeV changes from $13.87$ fb to $11.85$ fb when varying $\mu_{R}$ from $\sqrt{s}/2$ to $\sqrt{s}$, consistent with experimental measurements in Eqs. (1) and (2). The cross section at $\sqrt{s}=10.52$ GeV is $14.36$ fb and $12.29$ fb when setting $\mu_{R}$ at $\sqrt{s}/2$ and $\sqrt{s}$, respectively, serving as a prediction for the upcoming Belle II measurement. Finally, the total cross sections at $\sqrt{s}=91.2$ GeV and 250 GeV are respectively $8.48\times 10^{-7}$ fb and $2.84\times 10^{-10}$ fb, which could be of interest for future colliders.

Fig. 7 shows the total cross section for $J/\psi+\eta_{c}$ production, at the range of $\sqrt{s}=[6.5,12]$ GeV. The red, blue, and green bands stand for the LO, NLO, and NNLO cross sections, respectively. The width of each band represents the theoretical uncertainty from the renormalization scale variation, evaluated between $\mu_{R}=\sqrt{s}/2$ (upper edge) and $\mu_{R}=2\sqrt{s}$ (lower edge). A logarithmic scale is used for the vertical axis to clearly display the behavior near the physically interesting region around 10.6 GeV, given the cross section’s steep decline. The cross section peaks near 7 GeV, a value slightly above four times $m_{c}$. This comparison among the LO, NLO and NNLO results clearly demonstrates the significant contribution of the $\mathcal{O}(\alpha_{s}^{2})$ correction.

The total cross section over a broader energy range is presented on a logarithmic scale in Fig. 8. The cross section falls drastically with increasing energy. In this high-energy regime, higher-order corrections become increasingly significant due to the presence of large logarithmic terms.

VI.2 $\Upsilon+\eta_{b}$ production

Our analytical results and numerical asymptotic expansion retain explicit dependence on heavy quark mass, which makes them readily applicable to the process for the bottomonium counterpart. By substituting the relevant bottom quark parameters, as listed in Tab. 1, we obtain predictions for the cross section of $\Upsilon+\eta_{b}$ production.

From the numerical results shown in Tab. 3, we observe that the absolute magnitude of the cross section at $\sqrt{s}=22.5$ GeV is comparable to that of $J/\psi+\eta_{c}$. It can be seen from Tab. 2 and 3, the relative error of $J/\psi+\eta_{c}$ are generally larger than that of $\Upsilon+\eta_{b}$. For example, at $\sqrt{s}=$91.2 GeV, the relative uncertainty is approximately 28% for $J/\psi+\eta_{c}$ and 15% for $\Upsilon+\eta_{b}$. Additionally, for both processes, the relative error from the renormalization scale and the quark pole mass tend to decrease at higher energies.

Fig. 9 shows the cross section of $\Upsilon+\eta_{b}$ in the energy region from 20 to 150 GeV. Both axes are plotted on a logarithmic scale to enhance visibility of the high-energy region. Each band represents the variation associated with the renormalization scale $\mu_{R}$, where the upper and lower edges correspond to $\mu_{R}=\sqrt{s}/2$ and $2\sqrt{s}$, respectively. The cross section peaks around 22.5 GeV, where the cross section is

[TABLE]

As shown in Fig. 9, when $\sqrt{s}<30$ GeV, the NNLO corrections (the $\mathcal{O}(\alpha_{s}^{2})$ terms) contribute a relatively small fraction to the total cross section; this is visually manifested by the NNLO band (in green) being almost entirely covered by the NLO band (in blue). This reflects significantly better perturbative convergence compared to the $J/\psi+\eta_{c}$ case. Moreover, the uncertainty originating from the renormalization scale $\mu_{R}$ is noticeably reduced compared to NLO. However, when $\sqrt{s}>30$ GeV, NNLO corrections become more significant, which should be caused by large logarithms. It is worth noting that in the region where $\sqrt{s}\gg 4m_{Q}$, the cross section of $\Upsilon+\eta_{b}$ is much larger than that of $J/\psi+\eta_{c}$ due to the relatively larger long-distance matrix element.

VI.3 Convergence of the expansion in $r$

In this subsection, we analyze the convergence behavior of the two-loop amplitude $\mathcal{A}^{\mu(2)}$ within the expansion variable $r$ across different energy regimes. This is done by examining the contribution of the interference term

[TABLE]

For most phenomenological evaluations in this work, the asymptotic expansion is truncated at $\mathcal{O}(r^{35})$.

Fig. 10 displays the precision of partial cross section $\sigma_{20}$ with the asymptotic expansion up to $r^{i}$ with $i=2,10,12,14,16$. It can be seen that keeping two more powers of expansion brings one more efficient digit of precision. To quantify convergence, we define two functions describing the behavior as the expansion order $i$ increases:

[TABLE]

Here, $\sigma_{20}(i)$ is the partial cross section evaluated with the two-loop amplitude truncated at $\mathcal{O}(r^{i})$. The function $d(r,i)$ measures the actual precision achieved by the asymptotic expansion. The function $p(r,i)$ represents the predicted precision based on the convergence radius $1/16=0.0625$ derived from differential equation theory; it exhibits a linear dependence on $i$.

We analyze the asymptotic expansion for $r=0.01,0.02,0.03$ and $0.04$, corresponding approximately to center-of-mass energies of $15.0,10.61,8.66$ and $7.5$ GeV for $J/\psi+\eta_{c}$ production. Fig. 10 shows excellent agreement between $d(r,i)$ and $p(r,i)$. We find that achieving $6$-digit precision requires truncating the expansion at orders higher than $r^{8}$, $r^{12}$, $r^{18}$, and $r^{27}$ for $r=0.01,0.02,0.03$ and $0.04$, respectively. This confirms the predicted convergence behavior in the relevant energy range.

VI.4 Higher-order $\alpha_{s}$ contribution

To analyze the relative importance of each order in $\alpha_{s}$ expansion, we decompose the cross section as:

[TABLE]

where the terms are of order $\mathcal{O}(\alpha_{s}^{0})$ (first term), $\mathcal{O}(\alpha_{s}^{1})$ (second term), $\mathcal{O}(\alpha_{s}^{2})$ (third and fourth terms), $\mathcal{O}(\alpha_{s}^{3})$ (fifth term), and $\mathcal{O}(\alpha_{s}^{4})$ (sixth terms), respectively.

We plot the contribution of each term for $J/\psi+\eta_{c}$ production in Fig. 11. Both $\mathcal{O}(\alpha_{s}^{2})$ components (two green dashed lines) are sizable compared to the LO contribution (red solid line). Furthermore, the $\mathcal{O}(\alpha_{s}^{3})$ interference term, $2\operatorname{Re}[\mathcal{A}^{(2)}\mathcal{A}^{(1)*}]$ (black dashed line), remains significant, suggesting that the full $\mathcal{O}(\alpha_{s}^{3})$ correction may yield non-negligible contributions.

Fig. 11 shows that near $\sqrt{s}=10.6$ GeV, the partial $\mathcal{O}(\alpha_{s}^{4})$ term (purple dashed line) is much smaller than all lower-order contributions. However, as shown in Fig. 11, at higher energies this $\mathcal{O}(\alpha_{s}^{4})$ term, while still smaller than the $\mathcal{O}(\alpha_{s}^{1})$, $\mathcal{O}(\alpha_{s}^{2})$, and $\mathcal{O}(\alpha_{s}^{3})$ terms, becomes substantially larger than the $\mathcal{O}(\alpha_{s}^{0})$ LO term. This behavior raises concerns about the reliability of fixed-order perturbative expansions and may call for resummation or alternative theoretical treatments.

VII Summary and outlook

In this work, we have computed the $e^{+}e^{-}\rightarrow J/\psi+\eta_{c}$ two-loop amplitudes within the framework of NRQCD. All divergences from different sources are systematically subtracted, yielding finite amplitudes which are presented both analytically and numerically as an expansion in the dimensionless variable $r=m_{c}^{2}/s$. We provide a detailed analysis of the resulting cross section across various center-of-mass energies. By replacing $m_{c}$ by $m_{b}$ and changing the active number of flavors, we also extend our semi-analytical result to $\Upsilon+\eta_{b}$ production.

Our predictions are in good agreement with existing Belle and BaBar measurements, and we outline prospects for future experimental tests at Belle II and proposed high-energy lepton colliders. The analytical structure of our result provides a valuable check for the resummation formalism of large logarithmic terms in the future. Our results are presented at the amplitude level, and thus are well-suited for studying interference with other channels.

Nevertheless, certain theoretical uncertainties remain unsatisfactory. These include the dependence on the renormalization scale $\mu_{R}$, the convergence behavior of the perturbative QCD expansion, and the extraction of long-distance matrix elements within the NRQCD formalism. These limitations highlight the need for continued theoretical and experimental efforts in precision quarkonium physics.

Note added: While finalizing this paper, we noticed a paper Li et al. (2025) presenting a similar calculation. In this paper, not all the coefficients of the logarithms are expressed in analytic form.

Acknowledgements.

The authors would like to thank P. Zhang and X. Liu for useful discussion. This work is supported in part by the National Natural Science Foundation of China (Grants No. 11875071, No. 11975029, No. 12075251, No. 12005117, No. 12321005, No. 12035007), the National Key Research and Development Program of China under Contracts No. 2020YFA0406400, the High-performance Computing Platform of Peking University. XG is supported by the United States Department of Energy, Contract DE-AC02-76SF00515. XC is supported by the Swiss National Science Foundation (SNF) under contract 200020_219367. The Feynman diagrams in this paper are drawn with the aid of FeynGame Bündgen et al. (2025).