Color-octet contributions for J/ψ inclusive production at B factories in soft gluon factorization

We have studied color-octet contributions for J/ψ inclusive production at B factories, i.e., e+e− → J/ψ(3PJ8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {P}_J^{\left[8\right]} $$\end{document},1S08\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {S}_0^{\left[8\right]} $$\end{document}) + Xnon−cc¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {X}_{\mathrm{non}-c\overline{c}} $$\end{document}, using the soft gluon factorization (SGF) approach, in which the J/ψ energy spectrum is expressed in a form of perturbatively calculable short-distance hard parts convoluted with one-dimensional soft gluon distributions (SGDs). The series of velocity corrections originated from kinematic effect can be naturally resummed in this approach. Short-distance hard parts have been calculated analytically to next-to-leading order in αs. Renormalization group equations for SGDs have been derived and solved, which resums Sudakov logarithms originated from soft gluon emissions. Our final result gives a upper bound for color-octet matrix elements consistent with that extracted from hadron colliders. This may relieve the well-known universality problem in the NRQCD factorization. As a comparison, we also analytically calculated short-distance hard parts in the NRQCD factorization, with Sudakov logarithms resummed by using soft collinear effective theory. The comparison shows that velocity corrections from kinematic effect, which have been resummed in SGF, are significant for phenomenological study. Furthermore, it is found that Sudakov logarithms originated from soft gluon emissions are very important, while it is not the case for Sudakov logarithms originated from jet function. Therefore, the partial Sudakov resummation in SGF has already captured the main physics.


Introduction
Although heavy quarkonium production has been widely studied in the nonrelativistic quantum chromodynamics (NRQCD) factorization [1], the underline mechanism is still under debate. The reason is that the NRQCD factorization can not provide a universal description of all quarkonium production data. In other words, long-distance matrix elements (LDMEs) in NRQCD are found to be not universal. It was argued in ref. [2] that the universality problem in NRQCD may be caused by the bad convergence of velocity expansion, which suffers from large high order relativistic corrections due to soft hadrons emission in the hadronization process. Resummation of these relativistic-correction terms will result in the so called soft gluon factorization (SGF) framework [2]. It was demonstrated in ref. [3] that the SGF is equivalent to the NRQCD factorization, but with a series of important relativistic corrections originated from kinematic effects resummed. As a result, the SGF approach should has a much better convergence in the velocity expansion, and thus may provide a reasonable description of heavy quarkonium production.
Besides exclusive processes [4], the SGF approach has been recently applied to calculate the fragmentation function of the gluon to a 3 S [8] 1 heavy quark-antiquark pair in ref. [5], which was expressed in a form of perturbative short-distance hard part convoluted with one-dimensional 3 S [8] 1 soft gluon distribution (SGD). With a NLO calculation of the short-distance hard part, the authors demonstrated that the SGF is valid at NLO level.

JHEP03(2022)202
The renormalization group equation of the 3 S [8] 1 SGD was derived and solved, which resummed Sudakov logarithms to all orders in perturbation theory. A comparison with gluon fragmentation function calculated in NRQCD factorization indicates that the SGF formula resums a series of velocity corrections in NRQCD which are important for phenomenological study.
(1.2) However, this upper bound is much smaller than the value of CO LDMEs extracted from hadron colliders [28][29][30][31][32], which challenges the universality of LDMEs. Because SGF has a better convergence in velocity expansion, in this paper we apply it to study the CO contributions of e + e − → J/ψ( 3 P [8] J , 1 S [8] 0 ) + X non-cc . The rest of the paper is organized as follows. In section 2, we give a short review of SGF formula. Especially, we introduce a new lower cutoff x min of momentum fraction and demonstrate that final result is insensitive to the value of x min . In section 3, we present perturbative calculation in SGF. We also discuss the renormalization group equation (RGE) of SGDs. In section 4, we present perturbative calculation in NRQCD, and use the soft collinear effective theory (SCET) [26,[33][34][35][36] to resum large Sudakov logarithms to the NLL accuracy.
A comparison of phenomenological results obtained in SGF and NRQCD factorization are presented in section 5 and a summary is given in section 6. In appendix A we list some analytical expressions in perturbative calculation. Finally, we discuss the x min dependence in appendix B.
We also use light-cone coordinates where a four-vector a can be expressed as where the subscript ⊥ denotes the perpendicular direction. Then scalar product of two four-vector a and b becomes We introduce a light-like vector l µ = (0, 1, 0 ⊥ ), so that a · l = a + . Then in the rest frame of J/ψ we have where P ψ is the momentum of J/ψ, M ψ is the mass of J/ψ and q is the half of the relative momentum of the heavy quark-antiquark pair in J/ψ, which is relate to the relative velocity v by v = |q|/m c , where m c is the charm quark mass. The SGF is equivalent to the NRQCD factorization but with a series of important relativistic corrections originated from kinematic effects resummed [3]. Beginning from the leading operator for a specific quantum number in NRQCD Lagrangian, e.g., ψ † χ (or ψ † σ i T a χ, ψ † ← → ∇ i χ, and so on), one can construct powers suppressed operators like ψ † ← → ∇ 2 χ, ∇ 2 (ψ † χ) or ψ † gE · σχ by inserting the relative derivative ← → ∇ 2 , the total derivative ∇ 2 , or the E and B fields. The first two kinds of insertions are originated from kinematic effects and they are chosen to be resummed. Using equations of motion, one can replace the relative derivatives ← → ∇ 2 and ← → ∇ 0 by total derivatives, which results in operators like ∇ n 1 0 ∇ 2n 2 (ψ † χ). Then using integration by parts, one can eliminate all operators except that with n 1 = n 2 = 0. The price to pay is the introduction of the relative momentum between physical quarkonium and the intermediate heavy quark-antiquark pair into shortdistance hard parts [3]. The final formula is the SGF proposed in ref. [2]. Note that, in the above derivation, one needs to introduce proper gluon fields to combine with spacetime derivatives to form gauge covariant derivatives.
In SGF, the production cross section of J/ψ in e + e − collisions can be expressed in following formula [2]:

JHEP03(2022)202
where x is the longitudinal momentum fraction defined as x = P + ψ /P + with P ψ denoting the momentum of J/ψ and P denoting the total momentum of the intermediate cc pair. Different from the original form, we have introduced a lower cutoff x min for the x integration. This is allowed because, if x is too small, the intermediate cc pair will emit hard gluons during the hadronization process, which effect is perturbatively calculable. In appendix B we demonstrate that the production cross section is insensitive to the value of x min . In other words, the x min dependence of the integration limit cancels with the x min dependence ofσ [nn ] . This is not surprised at all becauseσ [nn ] is determined by matching the both sides of the above factorization formula. Therefore, as a special choice, one can also set x min = 0 [2].
In eq. (2.1), µ f is the factorization scale, √ s is the center-of-mass energy of the e + e − system, and the quarkonium state created by a † ψ has standard relativistic normalization. σ [nn ] are perturbatively calculable short-distance hard parts that produce an intermediate J ,J z in the amplitude and the complex-conjugate of the amplitude, respectively, with c, c = 1 or 8 representing the color-singlet or color-octet state of the cc pair. F [nn ]→ψ are one-dimensional SGDs which are defined explicitly as [5] where Ψ denotes the Dirac field of heavy quark. The subscript "S" means that the field operators in the definition are the operators obtained in small momentum region. In additional, we define "S" to select only leading power terms in threshold expansion [5], that is the expansion in the limit (P − P ψ ) + = (1 − x)P + → 0.
In general the state n can be different from the state n , but for the case of producing a polarization summed J/ψ, there are constraints c = c , S = S , J = J , J z = J z and |L − L | = 0, 2, 4, · · · [2,37]. While in this work, we only consider the case n = n with n = 1 S [8] 0 or 3 P [8] J,λ , where are the most important color-octet contributions for J/ψ production at e + e − collision. The corresponding projection operators K n , which define the intermediate state n, are given by [2] The color operator is given by with gauge link Φ l (rb − )ā a defined along the l µ direction,

JHEP03(2022)202
where P denotes path ordering and A µ (x) is gluon field in the adjoint representation: . E J,λ are polarization tensors for 3 P [8] J,λ states, with the following summation rules, where d is the space-time dimension, the spin projection operator P µν is defined as For energy distribution, we can rewrite eq. (2.1) as The variables E ψ /x and M ψ /x in the hard part correspond to energy E cc and invariant mass M cc of the intermediate cc pair, respectively. In above factorization formula, the short-distance hard parts H [n] are determined by the matching procedure [2,5]. To this end, we replace the final-state J/ψ by an on-shell cc pair with certain quantum number n and momenta To project the final-state cc pair to the state n, we replace spinors of the cc by following projector

JHEP03(2022)202
where, for n = 1 S [8] 0 or 3 P [8] J,λ , the color operator and spin operator are given by The factor N 2 c − 1 is to average over color-octet states. We insert the perturbative expansions In ref. [5], an explicit NLO calculation shows that the above kind of velocity expansion has good convergence, and the lowest order in the expansion can give a very good approximation of the full result. Therefore, to simplify the perturbative calculation, we only consider the contribution of H [n] here, with

The SGDs and the evolution equations
To obtain short-distance hard parts H (0) [n] up to NLO, we need first to calculate the perturbative SGDs. One-loop correction for the perturbative SGDs can be derived by following the calculation details of color-octet 3 S 1 SGD given in [5]. As we keep only the leading order in the velocity expansion, we can expand m 2 c around M 2 ψ /4 before performing the loop integration and phase space integration. The calculation is straightforward and we get Here ∆ = 1 − 4m 2 c /M 2 ψ , the "plus" functions are defined in the standard way through where f (x) is a well-behaved regular function.
In SGF, the general form of RGEs for SGDs is given by [5] where the evolution kernel can also be expressed in the velocity expansion series

JHEP03(2022)202
Using above perturbative SGD results we can obtain evolution kernel K at LO in α s by matching both sides of the RGE, where For the convenience of discussion, here we give some details of solving the above RGEs. We take 1 S [8] 0 SGD as an example. Following the discussion in ref. [5], we rewrite the RGE as a function of the variable where ω 0 = (1 − y)/y. ω 0 (or ω 0 ) is actually the longitudinal momentum fraction of the emitted soft gluons, and it takes values between 0 and ∞. To solve the above equation, it is convenient to perform the Laplace transformation, which, together with its inverse, is given byf where the constant c is chosen to be larger than the real part of the rightmost singularity off (s). Performing a Laplace transform in eq. (3.7) we obtain a RGE in Laplace space

JHEP03(2022)202
Here we evolved the SGD from the initial scale µ r /ν to the scale µ r .F mod [ 1 S [8] 0 ](ν) is a model introduced to describe the nonperturbative effects at the initial scale µ r /ν. And the evolution function h 0 is given by with where The Landau singularity can be avoided by performing following replacement in eqs. (3.11) and (3.12) as discussed in refs. [38,39] where a is a parameter of order 1 but not smaller than 1. Such a replacement prevents χ 0 from entering the nonperturbative regime. And the nonperturbative effects would be compensated by the modelF mod [ 1 S [8] 0 ](ν). Finally, with the help of the numerical inverse Laplace transformation, we can transform eq. (3.10) back to momentum space, The evolution of 3 P [8] J,λ SGD can be solved similarly.

The perturbative differential cross sections
We will calculate the perturbative differential cross sections of e + e − annihilate to the cc pair in state 1 S [8] 0 and 3 P [8] . For P-wave channel, we only calculate the polarization-summed hard part defined by (3.16) We choose the e + e − center-of-mass frame to perform the calculation. At LO, the Feynman diagrams are given in figure 1. Replacing the spinors of cc in the amplitude with the projectors in eq. (2.12) and integrating over the solid angle of q in the cc rest frame, we can derive for the 1 S [8] 0 channel. Herē For convenience, we use S and P to represent the cc[ 1 S [8] 0 ] and cc[ 3 P [8] ] states, respectively. For the P-wave contribution, we list the result in appendix A.
We take the S-wave contribution as an example to describe the calculation at NLO. Some typical NLO Feynman diagrams are showed in figure 2. The differential cross section is given by where we have use Lorentz covariance and gauge invariance to relate the cross section in e + e − annihilation to the decay rate of a virtual photon. Therefore, M Born denotes the tree-level amplitude for γ * → cc( 1 S [8] 0 ) + g, M virtual denotes the one-loop amplitude for γ * → cc( 1 S [8] 0 ) + g and M real is the amplitude for real correction of γ * → cc( 1 S [8] 0 ) + X, where X can be two gluons, a light quark-antiquark pair or a ghost-ghost pair. means summation over the polarizations final states and initial state virtual photon. dΦ 2(3) means two-body (three-body) phase space which are given by where q γ is the momentum of γ * , k and k i (i = 1, 2) represent the momentum of emitted gluon, light quark or ghost, as showed in figure 2, and N s is the symmetry factor for final-state particles.
In the calculation, we use FeynArts [40] to generate Feynman diagrams and amplitudes, and then use in-house code to contract Lorentz indices and carry out traces of Dirac matrices. According to eqs. (2.16), (2.18), we expand m 2 c in the amplitudes around M 2 ψ /4 and neglect the terms of O(|q|) before performing loop integration and phase space integration. We use the reverse unitary technique [41][42][43] to transform the delta functions to propagator denominators, then phase space integration can be treated similar as loop integration. We use the integration-by-parts (IBP) method [44,45] and employ the package FIRE5 [46] to perform the reduction of loop integrals, which express cross sections as linear combinations of master integrals (MIs). To calculate these MIs, we can set up differential equations with respect to r [47,48], The boundary conditions at r = 1 can be calculated by using the method of region [49]. We then use the algorithm in [50] to transform eq. (3.24) to the -form [51]. By solving the system, we obtain the MIs which are expressed by Goncharov polylogarithms [52]. We further simplify the expressions by using the package PolyLogTools [53]. Combining the IBP coefficients and the MIs we obtain expressions of the virtual and real corrections.
The UV divergence in virtual correction can be removed by the renormalization. We choose the renormalization constants Z 2 , Z m , Z 3 that correspond the charm quark field, the charm quark mass and the gloun field respectively in the on-mass-shell (OS) scheme, and JHEP03(2022)202 choose Z g that corresponds the QCD coupling in the minimal-subtractions (MS) scheme, Combining the virtual and real corrections, we get the differential cross section. The P-wave contribution can be calculated similarly. Expressions of S-wave and P-wave contributions are given by

JHEP03(2022)202
The expressions of R S , R P , R S and R P are listed in appendix A. Our analytical result for the S-wave differential cross section is equivalent to that in [54], and the analytical result for the P-wave channel is new.

Matching the short-distance hard parts
By inserting eqs.
One can find the short-distance hard part in P-wave channel is plagued with a singularity associated with the limit r = M 2 ψ /s → 1. Inserting the above hard parts into the SGF formula, M ψ will be replaced by M ψ /x, and thus such singularity appears at the threshold limit M 2 In the threshold limit, the emitted gluon at LO in figure 1 is very soft and its effect should be included in nonperturbative SGDs rather than in short-distance hard part. In fact, in the threshold region, 3 S [1] 1 channel can also contribute via e + e − → γ * → cc( 3 S [1] 1 ), followed by the hadronization process described by the SGD F [ 3 S [1] 1 ]→ψ . When 3 S 1 channel is included, the matching process to determine the 3 P [8] coefficient will have a contribution from the perturbative transition F [ 3 S [1] 1 ]→cc[ 3 P [8] ] , which will cancel the aforementioned threshold singularity.
However, in this paper we will be only interested in B-factories where r ≈ 0.1 is very small. Then, when x ∼ √ r, the nonperturbative quantities SGDs have significant perturbatively calculable effects, which can be avoided by introducing a reasonable x min . By choosing x min > √ r, both the contribution from e + e − → γ * → cc( 3 S [1] 1 ) and the threshold singularity of P-wave channel will disappear. This is consistent with the fact that SGF formula is insensitive to small x min .
At NLO, we have where T x = T | M ψ →M ψ /x . By integrating over x, we then obtain with θ(t) = 0 (1) if t is false (true), and We can find the above result still contains large logarithms at the limit z → 1 + r. These remained large logarithms originate from the collinear radiations recoil against the cc pair in the threshold region, which are not factorized in pure SGF. They can be resummed by introduce a jet function, as we will explain later.
Similarly, for the P-wave contribution we have

JHEP03(2022)202
We can find in the last line of eq. (3.33) that the singularity at z → 2 √ r is avoided due to the introduction of cut off x min > √ r. As a result, the cross section σ J/ψ is x min dependent. However, as demonstrated in appendix B, for small and moderate x min , the x min dependence of σ J/ψ is either α s suppressed or Λ QCD /M ψ suppressed. Thus the cross section σ J/ψ in SGF is insensitive to the choice of x min .

NRQCD factorization
As a comparison, we also present here the results of NRQCD factorization. The CO contribution for inclusive J/ψ production can be factorized as where the factor 1/9 is to average over polarizations of intermediate states, and O J/ψ ( 3 P [8] ) is the polarization-summed LDME defined by which has following relation at the lowest order in velocity approximation, Based on perturbative differential cross sections calculated in section 3.2, we can easily obtain the SDCs for NRQCD factorization and Using the same multi-loop calculation techniques, we can also obtain integrated cross sections. The corresponding SDCs are given bŷ where G S and G P are given in appendix A. Our results can reproduce that in ref. [24] if we choose the same parameters therein. Clearly, at the endpoint region, z → 1 + r (i.e. z → 1), the fixed-order results of the SDCs in eq. (4.4) suffer from large threshold logarithms. In order not to spoil the convergence of perturbative expansion, these threshold logarithms have to be resummed to all orders. Such resummation of the threshold logarithms to LO+NLL accuracy has been studied in ref. [26] within the SCET framework. According to [26,54], at the endpoint region we have following factorization formula where a phase space factor P [r , z ] is introduced [26] with P [r , 1] = 1. Fixed-order expression of the SDCs at the endpoint can be read from eq. (4.4), (4.9)

JHEP03(2022)202
The shape function S [n] (ξ) are defined in terms of ultrasoft fields that carry O(Λ QCD ) momentum where the ultrasoft covariant derivative is given as r), and ψ and χ denote the Pauli spinor fields in NRQCD that annihilates a heavy quark and creates a heavy antiquark, respectively. The jet function describes the collinear radiations recoil against the J/ψ in the threshold region, which is defined as (4.11) where the superscript (0) denotes the bare field, B µ ⊥ is the collinear gauge invariant effective field [26], and the factor s(1 + r) is chosen to provide a convenient normalization for the process considered here. According to [55,56], the fixed-order expression of the shape function and jet function are given by where M 2 J ≡ s(1 + r). Similar to SGF, to perform the resummation of the factorization formula eq. (4.8), we transform it to Laplace spacẽ

JHEP03(2022)202
where we introduce ω = 1−z and ω = 1−ξ, with ω denoting the longitudinal momentum fraction of the emitted ultrasoft gluons in the shape function. Different from the momentum fraction ω 0 = (1 − x)/x in SGD, which takes values between 0 and ∞, ω takes values between 0 and 1. Actually, the exact form of ω should be ω = (1 − ξ)/ξ, which also runs from 0 to ∞. However, in NRQCD+SCET approach one only consider the resummation in the endpoint region (ξ → 1), where one has ω = (1 − ξ)/ξ ∼ 1 − ξ, i.e. the evolution of the terms at higher order in 1 − ξ is not included. While in the evolution of SGD, these contributions are included. In eq. (4.13),S [n] (ν, µ) andJ (ν, µ) are the Laplace transformation of shape function and jet function, The componentsS [n] ,J and H [n] satisfy the renormalization group equations where anomalous dimensions obey relations Γ S cusp = −Γ J cusp = Γ H cusp and γ S + γ J + γ H = 0, and they are universal series in α s , Up to NLL accuracy, the needed coefficients are [26,57]

JHEP03(2022)202
We choose characteristic scales for H [n] , S [n] and J as follow then by solving the RGEs in eq. (4.15), we can evolve these functions from their characteristic scales to the scale µ r . Therefore, we get the resumed result where (4.20) and the evolution function is given by with Similar to the SGD, to deal with the Landau singularity, we make following replacement in eq. (4.21) (4.24) with branch points ν L 1 and ν L 2 given by . (4.25)

JHEP03(2022)202
The resummed cross section is modified as On the other hand, from eq. (4.9) we can derive the fixed-order expression ofσ[n](ν) as Comparing with eq. (4.26), we can find the single and double logarithms of ν in eq. (4.27) have been resummed. It should be noted that our choice of the characteristic scales in eq. (4.18) is different from that in refs. [26,54]. Their scales choices not only resum the logarithms of ν, but also try to resum logarithms of ln(s/m 2 c ). However, terms of ln(s/m 2 c ) can only be resummed by using the double-parton fragmentation framework [58][59][60][61], which is beyond the scope of this paper.
Using the numerical Laplace inverse transformation, we then obtain the resummed cross section in momentum space as where the subtraction of the last term is to avoid double counting between the endpoint and the full z range. As discussed before, in the pure SGF, only logarithms coming from soft gluon emission are resummed by using the RGEs of SGDs. And the evolution kernels for RGEs are known to O(α s ). For a comparison, here we also give two partly resummed results in NRQCD+SCET framework, which we denote as S-resum cross section and J -resum cross section. In deriving the S-resum cross section, the evolution of jet function is closed by choosing initial scales as Similarly, for J -resum cross section, the evolution of shape function is closed by following initial scales choice (4.31)

JHEP03(2022)202
In addition, in these two cases we only include the O(α s ) terms of the anomalous dimensions in eq. (4.16). Then these two partly resummed cross sections read

Phenomenology
In this section, we present phenomenological analysis based on our calculations in previous sections. The center-of-mass energy is chosen as √ s = 10.6 GeV for B factories. The QED

JHEP03(2022)202
coupling constant is set as α = 1/137. We determine the value of the coupling constant α s (µ r ) by adopting the two-loop RGE formula and setting Λ MS = 388 MeV, and we set the renormalization scale µ r = √ s/2 = 5.3 GeV. We take the J/ψ mass M ψ = 3.1 GeV and the heavy quark mass m c = 1.5 GeV. For nonperturbative modelF mod [n](ν) andS mod [n] (ν), we adopt the model function used in [26,54,55,62], which in Laplace space are given bỹ Following ref. [26], we assume that the parameters A, B andΛ in the model functions are the same for both of the S-wave and P-wave contributions. Furthermore, we also assume that these parameters are the same for the models of SGDs and shape functions. N ψ [ 1 S [8] Then after integrating out x, both the differential cross sections in SGF approach and in NRQCD approach can be written in the form And for the integrated cross sections we have In figure 3 we show the coefficients of the differential cross section in NRQCD factorization approach with different resummation methods. We can see that the resummed differential cross sections are much softer than the reult in NRQCD factorization. Moreover, in the NLO+NLL case the unphysical enhancement near the endpoint is cured by taking the resummation and nonperturbative shape function model into account. But the existence of the next-to-leading power logarithms, like ln(1 − z ), which are not resummed, still drives the differential cross section divergent in the endpoint region. We also find the S-resum cross section is close to the result in NLO+NLL method, while the J -resum cross section is deviate from the NLO+NLL result seriously. This phenomenon indicates that, as a good approximation, it makes sense to neglect the evolution of the jet function, as was done in pure SGF.
In figure 4 we show the differential cross sections in SGF comparing with that in NRQCD factorization. Due to the introduction of x min , the z distribution in SGF approach is piecewise function. Therefore, here we calculate the z distributions in 20 bins, and plot the average value for each bin. We find the peak in NRQCD+NLL results is  on the left of that in SGF. We also find that, comparing to the NLO+NLL results, the SGF results are significantly suppressed at moderate z, especially for the P-wave case.
There are many origins for these differences. First, in solving the RGEs of SGDs, we have included the effects at higher order in 1 − x, which will shift the peak to the right. While in the NRQCD+NLL method these effects have been ignored. Second, as shown in eqs. (3.29), (3.31) and (3.33), in SGF approach the hard parts in S-wave and P-wave are proportional to an overall factor x and x 3 (due to the factor 1/M ψ and 1/M 3 ψ ) respectively, which originated from velocity corrections resummation. Such factors suppress the differential cross sections at moderate z, especially for the P-wave. This implies that, even though large logarithms are resummed in NRQCD+SCET approach, there are still significant velocity corrections. Resumming these velocity corrections is the main purpose of the SGF. Finally, the SGF results included partial contributions at order O(α 3 s ), which come from the convolution of NLO short-distance hard parts with the resummed SGDs.
In table 1 we list the coefficients of the integrated cross section in different methods. We find the S-resum reults are close to the NLO+NLL reults. Besides, both the NLO+NLL results and SGF results are much smaller than the NRQCD reults due to resummation effects. Comparing to the NLO+NLL results, the SGF results are further suppressed.
In figure 5   moderate z as far as x min is not too large, e.g., x min < 0.55. In table 2, x min dependence of coefficients of cross section is shown. We find differences are smaller than 2%, which confirms our argument in appendix B that the cross section σ J/ψ in SGF is not sensitive to the choice of x min . Following ref. [24], we define a linear combination of CO LDMEs as: where "X" denotes factorization approach. If we set the CS contribution to be zero, and use the CO contribution to saturate the observed production cross section [9] σ[e + e − → J/ψ + X non-cc ] = 0.43 ± 0.13 pb,

JHEP03(2022)202
On the other hand, the value of CO LDME M k extracted from hadron colliders reads [28] M NRQCD,pp 3.9 = (7.4 ± 1.9) × 10 −2 GeV 3 , (5.8) which is about 3 times larger than the upper bound of M NRQCD 3.9 . Such a large discrepancy challenges the universality of NRQCD LDMEs. The SGF approach significantly reduces this discrepancy, which provides a hope to solve the universality problem. To this end, we must also describe J/ψ production in hadron colliders using SGF, which will be studied in future works.

Summary
In summary, in this paper we studied the J/ψ production via color-octet channel, e + e − → J/ψ( 3 P [8] J , 1 S [8] 0 ) + X non-cc , in SGF approach, in which a series of important velocity corrections can be resummed naturally. The corresponding J/ψ energy spectrum is expressed as a convolution of perturbatively calculable short-distance hard parts with one-dimensional color-octet SGDs, as shown in eq. (2.8). We introduced a cutoff x min for the longitudinal momentum fraction of the emitted soft gluons to prevent the gluons to be hard. We demonstrated both analytically and numerically that the x min dependence is suppressed and can be ignored in the sense of perturbation theory.
We calculated short-distance hard parts analytically up to NLO in α s . We derived and solved RGEs of SGDs, which resums Sudakov logarithms originated from soft gluons emission. Then by adopting a simple model for SGDs at an initial scale, our results are well-behaved near the kinematic endpoint, and they have the same shape as experimental data for the energy spectrum. By ignoring color-singlet contribution and using color-octet contributions to saturate the observed production cross section σ[e + e − → J/ψ + X non-cc ], we get a upper bound for the color-octet matrix element M X k , which seems to be consistent with the value extracted from hadron colliders.
However, there are two questions needing to be understood before we can convince ourselves that SGF has provided a reasonable description of J/ψ production. The first question is how important of the Sudakov logarithms originated from jet functions which have not been resummed in SGF. If this effect is very important, we need to do further job to resum them. The second question is how important of the velocity-correction terms that have been resummed in SGF. If this effect is not important, we do not really need the SGF approach, but simply using the NRQCD factorization combining with resummation methods to deal with Sudakov logarithms.
To understand the above two questions, we also calculated the same quantity in NRQCD factorization and use NRQCD+SCET approach to resum all encountered Sudakov logarithms. It was found that the full NRQCD+SCET result is very close to the result in which Sudakov logarithms originated from jet functions are not resummed. But if Sudakov logarithms originated from soft gluons emission are not resummed, one will get a result with large deviation from the full result. This answered the first question that Sudakov logarithms originated from jet functions are not that important. By comparing the

JHEP03(2022)202
NRQCD+SCET result with SGF result, we still find large difference, which implies that the resummation of velocity-correction terms is significant for phenomenological study.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.