Pseudoscalar Quarkonium+gamma Production at NLL+NLO accuracy

We consider the exclusive pseudoscalar heavy-quarkonium (eta_{b,c}) production in association with a photon at future lepton colliders where the collider energies of O(10^2) GeV are far greater than the quarkonium mass. At these energies, the logarithm of mass to collision energy becomes increasingly large hence its resummation becomes particularly important. By making use of the light-cone-distribution factorization formula, we resum the logarithms up to next-to-leading-logarithmic accuracy (NLL) that corresponds to order-alpha_s accuracy. We combine the resummed result with a known fixed-order result at next-to-leading order (NLO) such that both resummed-logarithmic terms and non-logarithmic terms are included at the same order in alpha_s. This allowed us to provide reliable predictions at accuracies of order alpha_s ranging from relatively low energies near quarkonium mass to the collider energies of O(10^2) GeV. We also include the leading relativistic corrections resummed at leading-logarithmic accuracy. Our prediction at the Belle energy is comparable with fixed-order predictions in literatures while it shows a large deviation from a recent Belle's upper limit by about 4 sigma. Finally, we make predictions for the energies of future Z- and Higgs factories.}


Introduction
Rigorous quantitative understanding of the heavy-quarkonium production at high-energy colliders [1] is a key probe not only to the features of quantum chromodynamics (QCD) but also to fundamental phenomena such as quark-gluon plasma (QGP) in heavy-ion collisions [2] and heavy-quark Yukawa couplings to Higgs [3,4]. An effective field-theoretic framework called the nonrelativistic QCD (NRQCD) [5] can be employed to predict quarkonium productions at high-energy colliders in a systematic way. NRQCD describes the dynamics inside a quarkonium at the energy scale m Q v 2 , where m Q is the mass of the heavy quark Q and v is the relative velocity of the Q andQ in the bound state. NRQCD is blind to the short-distance dynamics at higher energy scales of order m Q and the corresponding short-distance coefficients can be determined by matching to the full theory, QCD, which is known to be correct in all accessible energy scales. As a result, the production cross sections or decay rates involving heavy quarkonia can be expressed as linear combinations of NRQCD long-distance matrix elements (LDME) with the short-distance coefficients.
Exclusive processes such as associated photon production or double-quarkonium production at the lepton colliders like B factories and BES have been extensively studied in the framework of NRQCD. Future lepton colliders such as ILC [6], CEPC [7], and FCC-ee [8] offer opportunities to test our understanding of their productions at higher energies of O(10 2 ) GeV. In a collision at such a large center-of-momentum (CM) energy √ s, the cross section of a quarkonium has an uncomfortably strong dependence on the large logarithm of the ratio A straightforward extrapolation of the prediction for a lower-energy process of √ s 10 GeV to higher-energy processes listed above may result in a failure of predictive power. Thus the accuracy of a prediction can be reasonably controlled only after resumming the large logarithmic contributions in a proper way because such a logarithm cannot be suppressed by the strong coupling constant: α s ln r ∼ O (1).
The resummation of such a logarithm can be made by employing the light-cone (LC) approach [9,10] or, equivalently, the soft-collinear effective theory (SCET) [11]. In SCET, the scattering amplitude or current-current correlator is factorized into the following factors: the hardscattering kernel involving scales of √ s, the light-cone distribution amplitude (LCDA) that represents the collinear part, and the decay constant that involves the interactions of scales m Q . By solving the renormalization-group (RG) equation for collinear part or the hard part, one can resum the logarithms ln r.
In general, the collinear part describing a light meson such as a pion, ρ, or η is nonperturbative and one usually introduces an LCDA with a few model parameters. However, in the case of heavy quarkonium, the collinear parts can further be factorized into perturbative short-distance coefficients at the scale m Q and nonperturbative long-distance matrix elements at the scale m Q v 2 in the framework of NRQCD [12][13][14][15]. Therefore, it is worth to revisit and to update predictions by including the resummation of the large logarithms at energies of future lepton colliders. We express our formula in such a way that our expression reproduces the fixed-order results at low energies ∼ m Q , and at higher energies it resums large logarithms so that the same expression can be used for both the Belle and future high-energy experiments.
We consider the charge conjugate even (C = +1) processes with S-wave pseudoscalar quarkonium such as η c + γ and η b + γ. In a fixed-order perturbation theory this process was first computed in [16] at leading order (LO) and its next-to-leading order (NLO) correction was computed analytically in [17,18] and numerically in [19]. Up to date the α 2 s correction is available [20,21]. The relativistic correction of the order v 2 was first considered in [22] and α s v 2 correction was also obtained in [23]. The virtual Z contribution was computed up to α s correction in [24,25].
The leading-logarithmic (LL) accuracy resumming α n s log n r terms was first achieved in [12]. The quarkonium LCDAs at NLO were obtained by matching QCD onto NRQCD in [14,15]. In Higgs or Z boson decay into J/ψ + γ processes [26][27][28][29], the next-to-leadinglogarithmic (NLL) accuracy resumming α n+1 s log n r was achieved and the Abel-Padé method which enables to handle divergences appearing in computing the relativistic correction to the rates, was developed as well. Using the method, we make the prediction for η c,b + γ in lepton colliders at NLL+NLO plus the leading v 2 -correction accuracy.
The rest of the paper is organized as follows. In Sec. 2 we explain the theoretical formula to achieve NLL+NLO accuracy and provide all the ingredients for that order. Section 3 presents numerical results for the cross section and for the Z-boson decay rate into this process and Sec. 4 compares our result at the Belle energy to the previous results and to Belle's recent limit [30]. We finally summarize in Sec. 5.

Theoretical formula
The LC approach allows us to capture and to resum all logarithmic terms (singular), while non-logarithmic terms (nonsingular) can be computed by NRQCD fixed-order perturbation theory. We can express our full cross section as a sum of singular and nonsingular parts as in [31,32]. σ(r; µ, µ 0 , µ ns ) = σ sing (r; µ, µ 0 ) + σ ns (r, µ ns ) , where sing and ns in the superscripts and subscripts denote singular and nonsingular, respectively. In the singular part the scattering amplitude is factorized into a hard scattering kernel, an LCDA, and an NRQCD LDME. Each of them depends on a relevant energy scale such as µ ∼ √ s, or µ 0 ∼ m Q . 1 The renormalization group (RG) evolution between those scales enables us to resum large logarithms appearing in the fixed-order cross section and details of the evolution will be presented in coming subsections. On the other hand, if we turn off the resummation by setting all the scales being the same, it reduces to the singular part of the fixed-order cross section: σ fixed-sing (r; µ) = σ sing (r; µ, µ). The singular part of fixed-order cross section is given by where the coefficients are The born cross section σ 0 is given by where α is the fine structure constant, e Q is the fractional charge of a heavy quark Q, m P is the quarkonium mass, m Q is the heavy quark mass, O 1 P = 0|χ † ψ|P P |ψ † χ|0 is the non-relativistically normalized NRQCD LDME of the production of S-wave pseudoscalar quarkonium P , andẽ 2 Q is given in Eq. (2.18) which includes the effect of both of the virtual photon and Z boson propagators. To make our paper self-contained we also copy the fixedorder cross section in [17] into App. A.
The nonsingular part is defined by subtracting the fixed-order singular part from the fixed-order cross section as where the coefficients are c ns = (c fixed − c sing )/(−r) ≈ −14.9 − 4.8 log r + log 2 r + O(r) and Note that we pull out a prefactor −r in Eq. (2.5) to imply the relative suppression of nonsingular part in small r region but this correction is still important at the Belle energy. The nonsingular part from the Z-boson contribution is about 4m 2 Q /m 2 Z and remains small near the resonance and we omit them in this paper. Now the full cross section reproduces the ordinary fixed-order result when we turn off the RG evolution: σ fixed (r; µ) = σ(r; µ, µ, µ). Therefore, at the energies where log r ∼ O(1) hence µ ∼ µ 0 , the full cross section is consistent with the fixed-order results and at higher energies where | log r| 1 and µ µ 0 , the resummation implemented in full cross section becomes effective. Therefore the formula in Eq. (2.1) gives correct results at wide range of energies of current and future colliders. For the precise predictions for various energies including current and future colliders, both of the resummation of large logarithms and the fixed-order computation should be improved equivalently. Now let us discuss the amplitude for the singular part. In lepton collisions an initial lepton pair annihilates into virtual gauge bosons such as γ * or, Z * , then a pair of quark and anti-quark produced from the bosons turns into a bounded quarkonium state by emitting a photon. The scattering amplitude for a pseudoscalar quarkonium P plus a final photon can be written as where the index I = γ * , Z * represents the virtual bosons. The current has the vector and axial vector components: J µ V =ψγ µ ψ and J µ A =ψγ µ γ 5 ψ. But the axial current does not contribute due to the opposite charge conjugation and we only have J V in Eq. (2.6). The part L µ I contains a matrix element of initial lepton pair, a virtual boson propagator, and electroweak charges of the quarks. Its expression is where e is the electromagnetic coupling, e Q is the fractional electric charge of the heavy quark Q, √ s is the CM collision energy, θ W is the Weinberg angle, g e V = − 1 2 +2 sin 2 θ W and g e A = − 1 2 are the vector and axial charges of electron, and g Q V = T Q − 2e Q sin 2 θ W with T Q = ±1/2 is the vector charge of quark with a flavor Q = c, b, respectively.
The quark matrix element is extensively studied in the context of the meson form factor and the quark part in Eq. (2.6) can be expressed in the form factor style as where T is the time-ordered product and µν ⊥ = µνρσ p ρ k σ /(p · k) is the asymmetric tensor.

Light-cone distribution amplitude
The factor G P (µ) is the leading-twist result in LC factorization: where T P H is the hard-scattering kernel, φ P is the LCDA of a pseudoscalar quarkonium P , and f P is the decay constant. The scale µ is an arbitrary energy scale that separates the natural scales m Q and √ s which are the last arguments of each functions. Each function depends on the logarithm of their ratio: the natural scale to the scale µ. For simplicity, we omit the last arguments to the functions from now on.
The hard-scattering kernel T P H (x; µ) describes a production of quark and anti-quark pair at 1 S 0 state. The one-loop expression is given in Ref. [14,33] by where Note that ∆ = 0 for the naive dimensional regularization (NDR) scheme [34] and ∆ = 1 for the t'Hooft-Veltman (HV) scheme [35,36] for γ 5 regularization. 2 The scheme dependence in the hard kernel is cancelled by the same term with the opposite sign in the LCDA. The pseudoscalar LCDA is defined by a non-local matrix element of γ + γ 5 as where the plus components are γ + = γ 0 + γ 3 and p + = p 0 + p 3 . x ∈ [0, 1] is the collinear momentum fraction of a quark in the quarkonium and [z, 0] is the gauge link that is defined by where P stands for path ordering, g s = √ 4πα s is the strong coupling constant, A a is the gluon field with color index a, and T a is fundamental representation of SU(N c ).
The LCDA φ P describes the collinear-gluon exchange between quark and anti-quark pair and it is normalized to the unity upon the integration over x. This normalization defines the decay constant f P to be Eq. (2.12) at z = 0 and it describes a cc pair transition into a physical quarkonium. 3 In light mesons, the LCDA and the decay constant are nonperturbative and the former is modeled with a few parameters and the latter is determined by comparison to measurement. On the other hand in heavy quarkonium the LCDA and short-distance part of the decay constant are perturbatively calculable by matching QCD onto NRQCD amplitude. Their one-loop correction was obtained in [14] and the relativistic correction was obtained in [26,37]. We treat v 2 and α s corrections are of the same size and expand up to the same power. Then, the LCDA expanded up to the leading corrections in α s and v 2 is where , (2.15) and the + and ++ functions are defined in App. B. The leading α s and v 2 corrections to the decay constant f P are given by where Ψ P (0) is the wavefunction at the origin and is defined by 4 The relativistic correction agrees with the result in Ref. [37] with x 0 =x 0 = 1/2. To remove renormalon ambiguities, we replace the pole mass m Q in Eq. (2.16) with the MS mass m Q using the 1-loop relation [38]: and truncate higher-order contributions than our working precisions. The singular part of the cross section can be written in terms of G P (µ) in Eq. (2.9) as where for the virtual photonẽ 2 Q is e 2 Q , and for the virtual photon and Z boson it is is rather a fixed-order singular cross section in Eq. (2.2) because the functions Eqs. (2.10), (2.14), and (2.16) are in fixed-order form. We obtained the resummed singular part after the RG evolution and resummation, which is discussed next subsection.
We also give the expression for the decay rate of Z boson into a pseudoscalar quarkonium plus a photon in terms of G P (µ) as

RG equation and log resummation
The large logarithms in the cross section can be resummed by evolving each function in the factorization from its own natural scale, µ 0 ∼ m Q for the LCDA or µ ∼ √ s for the hardscattering kernel T P H to a common scaleμ, which can be chosen to be an arbitrary scale between µ 0 and µ because theμ dependence should be exactly cancelled when evolutions of all the functions are combined together. One of the simple and conventional choices is to setμ = µ then, the LCDA is just evolved from m Q to µ, while the hard-scattering kernel is treated as fixed-order function.
The LCDA evolution is governed by the RG equation called the Efremov-Radyushkin-Brodsky-Lepage (ERBL) equation [9,39]: where the ERBL kernel V (x, y; α s (µ)) for a pseudoscalar meson was extensively studied in the pion form factor. In the case of quarkonium the product f P φ P is factorized into two parts LDME and short-distance coefficients as in Eqs. (2.16) and (2.15) and the RG equation Eq. (2.20) can be expressed into two set of RG equations: one for LDMEs and the other for the coefficients. The formal equation is evolved from LDME's natural scale m Q v 2 to µ and the latter is from µ 0 to µ. This way would better fit to the philosophy of scale separation in effective field theory framework. However, LDME scale is nonperturbative and evolution from the scale would not work. Conventionally the LDMEs are determined at the scale µ 0 rather m Q v 2 . Then, we simply run both parts from µ 0 to µ by using Eq. (2.20).
The LL and NLL accuracies are achieved by solving the ERBL equation with one-and two-loop kernels respectively. The kernel is known up to two loops in the NDR scheme and we use the NDR results to achieve NLL accuracy.
where the one-loop coefficient is given by 22) and the two-loop expression can be found in Refs. [40][41][42][43]. The eigenfunction of the one-loop kernel where G n is the product of the Gegenbauer polynomial C (3/2) n and its weight [44]: n is the LO anomalous dimension and here we follow the convention of Ref. [45] The NLO anomalous dimension γ (1) n is defined in the same way from the 2-loop kernel V (1) (x, y). The solution of the ERBL equation is expressed as a series sum, where the k-th Gegenbauer coefficient of LCDA at the scale µ 0 is for even n and zero for odd n. Following the convention of Ref. [45], the solutions of the scale evolution factor U nk (µ, µ 0 ) up to NLL accuracy are given by where the value of U nk at LL accuracy is nonzero only for n = k: [α s (µ)/α s (µ 0 )] γ (0) n /(2β 0 ) δ nk . At NLL it is nonzero when (n−k) is zero or, even and positive integer. β n is the beta function coefficients for (n + 1)-th order in α s and the explicit expressions of the two-loop anomalous dimensions γ (1) n and d nk are copied in App. C. Then, the RG evolved function G P (µ) is given by where M (i,j) are defined in terms of the RG evolved LCDA in Eq. (2.26) and is the coefficient in the expansion with Gegenbauer polynomials and it is non-vanishing only for even n because T P H is symmetric with respect to x = 1/2 while G n (x) is asymmetric for odd n. The LO coefficient for even n is simple We would like to note that the decay constant f P in Eq. . We also emphasize that the relativistic correction M (0,v 2 ) (µ) correctly resums logarithms proportional to v 2 α n s log n r by using the same RG evolution U nk and its expansion in α s is given by This agrees with the logarithmic term in α s v 2 correction Eq. (2.26) in [23]. Eventually we insert Eq. (2.30) into Eq. (2.17) and obtain the singular part and full cross section in Eq. (2.1). In practice of computation there is an option of truncating higher-order terms irrelevant at our NLL accuracy and we make following truncations. In Eq. (2.30), the first term M (0,0) is computed using NLL expression of U nk while the other M (i,j) terms are computed using LL expression. In the absolute square |G P (µ)| 2 , we also drop higher-order terms proportional to α 2 s or α s v 2 , which are obtained in the product of M (i,j) in Eq. (2.30) and f P in Eq. (2.16).
We also adopt the Abel-Padé method developed and used in [27][28][29] to achieve faster numerical convergence at NLL accuracy and to deal with divergences associated with the relativistic corrections in LCDA.

γ 5 -scheme dependence
As we can see from Eqs. (2.11b), (2.15), and (2.16), there are γ 5 scheme dependences in the hard kernel T H (x, µ), LCDA φ P (x, µ) and decay constant f P (µ), which are represented by the terms proportional ∆ = 0, 1 for NDR and HV schemes, respectively. It is easy to check that the ∆ dependences of the factor G P (µ) vanish at NLO without resummation or, RG evolution. However, it is not obvious whether the ∆ dependences vanish or not at NLL accuracy due to additional scheme dependence that may enter in two-loop anomalous dimension γ (1) n . Note that γ (0) n is ∆-independent and so is the LL resummation. Ref. [46] computed the n f -dependent part of two-loop evolution kernel V (1) (x, y) in both NDR and HV schemes and we can obtain the scheme dependence for full two-loop evolution kernel by combining Eqs. (5.24), (5.35), and (5.41a) and applying the relation in Eq. (5.40) in Ref. [46], which gives Again, polynomials G n (x) is the eigenfunction of ∆V (1) with the eigenvalue anomalous dimension where analytic expression of the anomalous dimension is given by . (2.37) At NLL, there are two types of ∆ dependences in the amplitude G P (µ). The one from NLL evolution factor U nk in Eq. (2.29) is proportional to the anomalous dimension in Eq. (2.37): Here G LL P,n (µ) = αs(µ) n are given in Eqs. (2.28) and (2.33). The other type of ∆ dependence is those from one-loop corrections f (1) n . The terms proportional to ∆ are given by non-logarithmic constant parts ∆f (1) Collecting three contributions above, we have (2.42) In the second equality we rearranged terms into two parts, the one proportional to NLO correction evaluated at the α s scale µ and the other proportional to the difference α s (µ 0 ) − α s (µ) which is a part of NLL resummation. The scheme independence at fixed-order NLO implies the cancellation of first part This is also confirmed by explicit computing ∆T (1) n and ∆φ (1) n which are zero for odd n and for even n. We would like to note an interesting relation between ∆T n we eliminated ∆f (1) n and ∆φ (1) n in favor of ∆T (1) n by using Eq. (2.43). Therefore, γ 5 scheme independence is valid at NLL accuracy.
At higher-order, we expect a similar pattern of cancellation between ∆ dependent terms: cancellation between fixed-order terms at the same α s scale as in Eq. (2.43) and cancellations between n-loop anomalous dimension from evolution factor and the constant terms of (n − 1)loop function as in Eq. (2.47).

Logarithmic structure
Here we discuss logarithmic structure and accuracy of the resummed amplitude. Even though this section explains quite well-known properties of resummation and does not contain anything new, it may be useful for those who are not familiar with resummation.
Let us first look at the fixed-order expansion of amplitude in Eq. (2.9). Its logarithmic structure can be schematically expressed as where L ≡ log(4m 2 Q /s). The largest logarithmic term at each α s order is α n s L n and then the next largest is α n s L n−1 . The functions T H and φ are first expanded with the Gegenbauer polynomials, then each coefficient is resummed as G nk P = f P T n U nk φ k , and summation over all n and k gives the resummed G P . The individual G nk P takes the following form where C(α s ) is the fixed-order expansion in α s and it does not depend on the logarithms, (2.50) The fixed-order coefficient C i is given by fixed-order function T n , φ k , f P and the coefficients C ij associated with anomalous dimensions are given by U nk in Eq. (2.29). Therefore, in general the coefficients C i and C ij differ for the different values of n, k. For example, C 0 is non-zero for the diagonal element where n = k but zero otherwise. We are implicit with those n, k dependence to make our discussion focused on the logarithmic structure. Similarly, we do not separately discuss about v 2 corrections in the coefficients C i and it follows the same conclusion.
In the fixed-order perturbation theory, the series in Eq. (2.48) are summed row-by-row, i.e., order-by-order in α s . On the other hand, in resummed perturbation theory, the series in the exponent of Eq. (2.49) are summed column-by-column, based on large-logarithmic power counting L ∼ 1/α s . In Eq. (2.49), the first column is of the order α n s L n ∼ 1 called the LL, the second column is α n s L n−1 ∼ α s called the NLL, and so on. It is clear that which fixed-order terms in Eq. (2.50) should be included: C 0 at LL and C 1 at NLL.
However, one may realize that the structure of evolution factor U nk in Eq. (2.29) is different from that of Eq. (2.49). For example, the non-exponent term contains logarithms: 2π β 0 log(µ/µ 0 ) + · · · . This is because, in Eq. (2.49) the second column is of O(α s ) hence those terms beyond LL can be expanded and moved down from the exponent: whereC(α s ) includes two prefactors (2.52) We have LL accuracy with first term in Eq. (2.52) and NLL with O(α s ) terms in the largelogarithmic power counting α n+1 s L n ∼ α s . This alternative way of arranging logarithms is equivalent to Eq. (2.49) up to higher-order corrections than working accuracy and is the formula we use in this paper.

Numerical results
In this section, we list input parameters for numerical calculations then, present our results for the final state η c,b + γ in e + e − collisions at various collision energies and in Z-boson decay. Those results include the resummation at NLL accuracy, the fixed-order correction at NLO, and the relativistic corrections of the order v 2 as we discussed in previous sections. The numerical results at LL and NLL+NLO are compared and their perturbative convergence is discussed. The uncertainty includes variations of σ, m c , m 2S − m 1S , Γ[η c → γγ]. And we assumed the size of the neglected higher-order corrections in α s and v 2 to be 30% times the central values of α s and v 2 , respectively. Major sources of uncertainty are the variation of σ and the assumed higher-order corrections. We like to pay a bit more attention to using those values in Eq. (3.1). In conventional predictions, we may use the same value for the predictions at LO, NLO, and higher-order accuracy. This way can correctly reproduce the input decay rate using its 1-loop expression, while with LO expression of the decay rate the result is systematically biased by the amount of 1-loop correction included in the matrix element. Of course it is not a problem when the size of 1-loop correction is small as in bottomonium. However in the case of η c we find the effect is as large as 40% due to large coupling constant α s (2m c ) ∼ 0.27. A similar bias exists in the cross sections and leads to overshooting in its predictions at LO. Eventually it would spoil the perturbative convergence due to a large change from LO to NLO and similarly from LL to NLL.

Input parameters and NRQCD matrix elements
We can avoid this systematic bias once we use the matrix element determined at the same order with working order, at which we make predictions. By doing this the (experimental) input decay rate is always reproduced at each order in α s . This can be done by replacing the NRQCD matrix element with the experimental value of decay rate multiplied by short distance coefficient: where the experimental value for η c is Γ exp = 5.0 ± 0.4 keV [51]. We insert Eq.
Note that the coefficient of α s C F /π reduces from −3 to (8 − π 2 )/4 ≈ −0.5 and the coefficient of v 2 changes from −2 to −2/3. In this way, we have better perturbative convergence between LL and NLL (LO and NLO).

Final results
Our numerical results for η c + γ cross sections and perturbative uncertainties at different accuracies are given in Fig. 1. Three accuracies LL, NLL combined with NLO non-singular part and leading v 2 correction (NLL+NLO), fixed-order NLO are compared. The bands on left and right panels are absolute and relative perturbative uncertainties. The cross section in figures is scaled by the LO cross section The values of scales we choose for LL and NLL+NLO are µ 0 = 2m c and µ = µ ns = √ s and for NLO we set all scales to be the same µ 0 = µ = µ ns = √ s. The perturbative uncertainties are estimated by varying µ, µ ns from its central value by a factor 2 up and down and by varying µ 0 by a factor of √ 2. The uncertainties are summed in quadratures as in [31]: δσ 2 µ 0 + δσ 2 µ + δσ 2 µns , where δσ µ i is the change of cross section by a scale variation of  Table 1. The perturbative uncertainty (width of the band) decreases by a factor of α s from LL to NLL+NLO and a reasonable overlapping between two bands in left panel implies a good perturbative convergence. With increasing CM energy, the deviation of NLO from NLL+NLO becomes more significant due to the large logarithms not taken into account at NLO and this clearly shows that the small perturbative uncertainty of NLO is not reliable at this high energies. In Fig. 2 we also show the results with a smaller value of µ 0 : µ 0 = m c instead of 2 m c . The NLL perturbative uncertainty at µ 0 = m c tends to be asymmetric and smaller than that for µ 0 = 2m c because the lower scale variation from m c by a factor of √ 2 moves µ 0 close to the Landau pole and the scale dependence near this region is not monotonic. In comparison to Fig. 1 at µ 0 = 2m c we observe relatively better perturbative convergence between LL and NLL although the other is still reasonable. For these reasons we take µ 0 = 2m c for our final results listed in Table 1.
In Fig. 3 we also show the bottomonium production cross sections and their percent perturbative uncertainties, which are smaller compared those for charmonium. The decay rate for η b is not available we use following value Γ[η b → γγ] = 0.512 +0.096 −0.094 keV and for relative velocity v 2 η b = −0.009 +0.003 −0.003 taken from [52]. The central values of scales are µ 0 = 2m b and µ = µ ns = √ s and their variations are done in same way as for the charmonium. Table 1 lists our final results for the cross sections at B-, Z-and Higgs-factory energies: √ s =10.58, 91.2, and 240 GeV and Z-boson decay branching fractions. The uncertainties in the table for charmonium (η c ) channel includes uncertainties of input decay rate Γ exp (±8%), relative velocity v 2 (±30%) as well as perturbative uncertainties (±3% or less) shown Fig. 1 and they are added in quadratures. For the bottomonium (η b ) the uncertainties are input decay rate (±20%), relative velocity (±30%), perturbation (±3% or less). The final uncertainties quoted in Table 1 Figure 4. Status of NRQCD predictions (points) and Belle's upper limit [30] (horizontal gray line) for e + e − → η c + γ. This work on right side for LL and for NLL+NLO with v 2 corrections is compared to previous predictions (from the left) at NLO [19], NLO with v 2 correction [17] and with v 2 resummation [22], and NNLO [21]. We note that uncertainties of NLO [19] and of NNLO [21] include only quark mass variations while the others include other sources of uncertainties. See the text for more details.

Comparison to various predictions and Belle's upper limit
In Fig. 4, we summarize the status of NRQCD predictions (points) in comparison with Belle's upper limit (90% credibility level) [30] (gray line) for σ(e + e − → η c + γ) at √ s = 10.58 GeV. Since the resummation effect is not substantial at this energy our results LL and NLL+NLO+v 2 on the right side of the plot should be comparable to LO and to NLO with v 2 corrections.

corrections.
For a fair comparison with previous predictions we need to point out several major differences of input parameters and their variations between different predictions. First, a small error bar of NLO [19] and invisibly small error of NNLO [21] only include a charmquark mass variation by 0.1 GeV and they should not be compared to full uncertainties of other predictions. Instead their central values can be compared with the others. Second, LO, NLO+v 2 [17], and NLO+v 2 resummation [22] use the LDME of [50], which should be updated with improved measurement of Γ[η c → γγ] as discussed around Eq. (3.1) and with the updated LDME, we expect decrease of the cross section by about 10 ∼ 20% and also reduction of their uncertainties, quantitative estimation of which requires more careful study and is beyond scope of this paper. On the other hand our results of LL, NLL+NLO+v 2 in Fig. 4 is lower in its value and smaller in uncertainty partially due to this update.
There are differences in scale choice and its variation. We use two scales µ = √ s for the hard-scattering kernel and µ 0 = 2m c for the LCDA and decay constant and they are varied by a factor of 2 for µ and √ 2 for µ 0 as discussed in previous section. While we use MS mass m c = 1.275 +0. 25 −0.35 , many of previous results use the pole mass m c = 1.4 ∼ 1.6 GeV. NLO [19] sets µ = 2m c and NLO+v 2 resummation [22] sets µ = m c , 2m c (result with 2m c is shown in Fig. 4 and the value for m c is similar), while NLO+v 2 [17] makes most conservative choice µ = √ s, 2m c , and m c , which leads to relatively larger uncertainty compared that of NLO+v 2 resummation. NNLO [21] chooses different values for renormalization scale µ r = √ s/2 and and the factorization scale µ Λ = 1.0 GeV. We do not include the result of [18] because the 1-loop coefficient is not consistent with other results [14,17].
Recently the Belle experiment analyzed S-wave (η c + γ) and P-wave (χ cJ + γ with J = 0, 1, 2) channels [30]. While P-wave cross section (J = 1) and upper limits (J = 0, 2) are consistent with the theoretical predictions [16][17][18][19], S-wave upper limit σ exp ηc+γ ≤ 21.1 fb at 90 % credibility level is in tension with our NLL+NLO prediction 32.7±2.8 fb by 4.1σ. This reminds us the puzzle in exclusive J/ψ + η c production [53][54][55][56], where a large discrepancy between theory and experiment was resolved by the combined effect of large K-factor, resummed relativistic corrections and careful determination of NRQCD matrix element [50,[57][58][59]. 7 However, in our case effect of the K-factor, a ratio of NLL+NLO including relativistic correction relative to LO, is less than 5% as shown in Fig. 1. It would be surprising if higher-order resummation or relativistic corrections is the resolution to this tension. Of course, more careful study on those corrections and other contributions from different topology can shed lights on the tension. Without correct understanding of this channel one may also cast a doubt on theoretical prediction for other exclusive processes such as radiative Higgs decay into quarkonium, a novel channel to probe the Yukawa coupling of charm quark [3,4,61,62]. In this aspect, resolving the tension would be one of important checkpoints. The Belle II experiment with upgraded luminosity is starting its physics program and in a few years it will release improved measurements and can clarify if this seemingly tension is to be or not.

Summary
We resum large logarithms of 4m 2 Q /s at NLL accuracy for exclusive production for η b,c + γ in high-energy lepton colliders by using light-cone factorization theorem and by using 2-loop evolution kernel known from the pion form factor. The leading relativistic correction is also included and logarithms in the correction is resummed at LL accuracy. The nonsingular part of order α s is obtained by subtracting the singular part from fixed-order results at NLO then, is added to resummed cross section. This makes our prediction of order α s accuracy valid in both resummation region (r 1) and fixed-order region (r ∼ O(1)) where r = 4m 2 Q /s such that the results with the same formalism in Eq. (2.1) can be compared to measurement at the Belle energy near 10 GeV and in future colliders such as ILC, CEPC, FCC-ee.
Our final state η c,b is the pseudoscalar, which involves an ambiguity in handling γ 5 in d dimension and the scheme dependency enters in individual parts such as hard kernel, LCDA, and the decay constant in factorized formula. We explicitly showed that how the γ 5 -scheme dependence vanishes in the resummed expression at NLL accuracy. In resummed expression, there is a part proportional to fixed-order singular result and its scheme independence is followed by that of the fixed-order cross section. In the other part of resummed expression, we observe that the scheme dependence of 2-loop anomalous dimension is matched to and cancelled against constant term of 1-loop hard-scattering kernel Eq. (2.47).
In numerical calculation in Sec. 3 we first rewrite the decay constant in terms of the experimental decay rate by eliminating NRQCD matrix element to avoid a systematic bias by unnecessary higher-order α s contribution that can be contained in NRQCD matrix element. By doing this all the input formula are computed at the same α s order to working accuracy and it is observed to show better perturbative convergence from LO to NLO and from LL to NLL. Our predictions for the cross sections and branching fraction are summarized in Table 1. The input decay rates Γ[η b,c → γγ] dominates over the others including perturbative uncertainty and uncertainty of our prediction reduces if the measurement of decay rate improves. In Sec. 4 we compare our prediction to previous predictions for the Belle experiment and discussed the difference input parameters and uncertainty estimate. We also find Belle's recent upper limit 21 fb is about 4 σ away from our prediction 33 ± 3 fb. We hope future Belle II analysis with better statistics coming-out in a few years and careful theoretical investigation on higher-order corrections may shed lights on this tension.