Pseudoscalar Quarkonium+gamma Production at NLL+NLO accuracy

TL;DR

This paper develops a theoretical framework to accurately predict pseudoscalar quarkonium plus photon production at high energies by resumming large logarithms and combining with fixed-order calculations, relevant for future collider experiments.

Contribution

It introduces a combined resummation and fixed-order approach for NLL+NLO accuracy in pseudoscalar quarkonium production, extending predictions across a wide energy range.

Findings

Predictions align with existing literature at Belle energies.

Significant deviation from recent Belle upper limits (~4 sigma).

Provides forecasts for future Z and Higgs factory energies.

Abstract

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.

Figures

Tables

Table 1. Table 1: Cross sections σ ​ ( e + ​ e − → η Q + γ ) 𝜎 → superscript 𝑒 superscript 𝑒 subscript 𝜂 𝑄 𝛾 \sigma(e^{+}e^{-}\to\eta_{Q}+\gamma) and branching fractions Br( Z → η Q + γ → 𝑍 subscript 𝜂 𝑄 𝛾 Z\to\eta_{Q}+\gamma ) for charmonium and for bottomonium with uncertainties including all input parameters as well as perturbative uncertainties.

	Cross section		Branching fraction
$\sqrt{s}$	${\eta }_c$	${\eta }_b$	${\eta }_c$	${\eta }_b$
10.58 GeV	$32.7\pm 2.8$ fb	-	$\left(7.42\pm 0.61\right)\times {10}^{-9}$	$\left(2.80\pm 0.53\right)\times {10}^{-8}$
$m_Z$	$0.449\pm 0.037$ fb	$1.66\pm 0.31$ fb
240 GeV	$0.189\pm 0.016$ ab	$0.0934\pm 0.0176$ ab

Equations

$$r\equiv \frac{4m_Q^2}{s}.$$

$${\alpha }_s\ln r\sim O(1).$$

$$\sigma (r;\mu ,{\mu }_0,{\mu }_{\text{ns}})={\sigma }^{\text{sing}}(r;\mu ,{\mu }_0)+{\sigma }^{\text{ns}}(r,{\mu }_{\text{ns}}),$$

$${\sigma }^{\text{fixed-sing}}(r;\mu )=\underset{r\to 0}{\lim }{\sigma }^{\text{fixed}}(r;\mu )={\sigma }_0\left[1+\frac{{\alpha }_s{C}_F}{4\pi }{c}^{\text{sing}}+⟨v^2⟩c_{v^2}^{\text{sing}}\right],$$

$$2m_P\left[1+⟨v^2⟩c_{v^2}^{\text{sing}}\right]=4m_Q\sqrt{1+⟨v^2⟩}\left[1+⟨v^2⟩c_{v^2}^{\text{sing}}\right]\approx {\left[2\sqrt{m_Q}\left(1-\frac{5}{12}⟨v^2⟩\right)\right]}^2,$$

$$c^{\text{sing}}$$

$$c_{v^2}^{\text{sing}}$$

$${\sigma }_0=\frac{16{\pi }^2{\alpha }^2(\sqrt{s})\alpha (0)e_Q^2{\tilde{e}}_Q^2{m}_P}{3s^2}\frac{⟨{\mathcal{O}}_1{⟩}_P}{m_Q^2},$$

$$⟨v^2⟩\equiv \frac{1}{m_Q^2}\frac{⟨P|{\psi }^{†}{\left(-\frac{i}{2}\overleftrightarrow{\mathbf{D}}\right)}^2\chi |0⟩}{⟨P|{\psi }^{†}\chi |0⟩},$$

$${\sigma }^{\text{ns}}(r;\mu )={\sigma }^{\text{fixed}}-{\sigma }^{\text{fixed-sing}}=-r{\sigma }_0\left[1+\frac{{\alpha }_s{C}_F}{4\pi }{c}^{\text{ns}}+⟨v^2⟩c_{v^2}^{\text{ns}}\right],$$

$$i\mathcal{M}(e^{+}{e}^{-}\to {\gamma }^{\ast }/Z^{\ast }\to P\gamma )=L_I^{\mu }⟨P(p)+\gamma (ϵ,k)|J_V^{\mu }|0⟩,$$

$$L_I^{\mu }=\{\begin{array}{ll}\frac{i{e}_Q{e}^2}{s}\bar{v}(\bar{P}){\gamma }^{\mu }u(P),& \text{for I=γ∗},\\ -\frac{i{g}_V^Q{e}^2/{\sin }^2(2{\theta }_W)}{s-m_Z^2}\bar{v}(\bar{P}){\gamma }^{\mu }(g_V^e-g_A^e{\gamma }_5)u(P),& \text{for I=Z∗},\end{array}$$

$$⟨P(p)+\gamma (ϵ,k)|J_V^{\mu }|0⟩$$

$$G_P(\mu )\equiv f_P(\mu ,m_Q){\int }_0^{1}dx{T}_H^P(x;\mu ,\sqrt{s}){\varphi }_P(x;\mu ,m_Q),$$

$$T_H^P(x,\mu )=T_H^{(0)}(x)+\frac{{\alpha }_s(\mu )}{4\pi }{T}_H^{(1)}(x,\mu )+O({\alpha }_s^2),$$

$$T_H^{(0)}(x)$$

$$T_H^{(1)}(x,\mu )$$

$$⟨P(p)|\bar{Q}(z){\gamma }^{+}{\gamma }_5[z,0]Q(0)|0⟩=p^{+}{f}_P{\int }_0^{1}dx{e}^{ip\cdot zx}{\varphi }_P(x,\mu ),$$

$$[z,0]=\mathcal{P}\text{exp}\left[i{g}_s{\int }_0^{z}dy{A}_a^{+}{T}^a\right],$$

$${\varphi }_P(x,\mu )={\varphi }^{(0)}(x)+\frac{{\alpha }_s(\mu )}{4\pi }{\varphi }^{(1)}+⟨v^2⟩{\varphi }^{(v^2)}+O({\alpha }_s^2,{\alpha }_s{v}^2,v^4),$$

$${\varphi }^{(0)}$$

$${\varphi }^{(1)}$$

$${\varphi }^{(v^2)}$$

$$f_P(\mu )=\frac{\sqrt{2N_c}\sqrt{2m_P}{\Psi }_P(0)}{2m_Q}\left[1-⟨v^2⟩+\frac{{\alpha }_s(\mu )C_F}{4\pi }(-6+4\Delta )+O({\alpha }_s^2,v^4,{\alpha }_s{v}^2)\right],$$

$${\sigma }^{\text{sing}}=\frac{2{\pi }^2{e}_Q^2{\tilde{e}}_Q^2{\alpha }^2(\sqrt{s})\alpha (0)}{3s^2}|G_P(\mu ){|}^2,$$

$${\tilde{e}}_Q^2=e_Q^2-2\frac{e_Q{g}_V^Q{g}_V^e}{{\sin }^2(2{\theta }_W)}\frac{1-r_Z}{(1-r_Z{)}^2+r_Z\frac{{\Gamma }_Z^2}{s}}+\frac{(g_V^Q{)}^2[(g_V^e{)}^2+(g_A^e{)}^2]}{{\sin }^4(2{\theta }_W)}\frac{1}{(1-r_Z{)}^2+r_Z\frac{{\Gamma }_Z^2}{s}},$$

$${\Gamma }^{\text{fixed-sing}}(r;\mu )=\frac{\pi \alpha (m_Z)\alpha (0)e_Q^2(g_V^Q{)}^2}{6m_Z{\sin }^{2}2{\theta }_W}|G_P(\mu ){|}^2.$$

$${\mu }^2\frac{\partial }{\partial {\mu }^2}\left[f_P(\mu ){\varphi }_P(x,\mu )\right]={\int }_0^{1}dyV(x,y;{\alpha }_s(\mu ))\left[f_P(\mu ){\varphi }_P(y,\mu )\right],$$

$$V(x,y;{\alpha }_s(\mu ))=\sum _{n=0}^{\infty }{\left(\frac{{\alpha }_s(\mu )}{4\pi }\right)}^{n+1}{V}^{(n)}(x,y),$$

$$V^{(0)}(x,y)=2C_F{\left[\frac{1-x}{1-y}\left(1+\frac{1}{x-y}\right)\theta (x-y)+\frac{x}{y}\left(1+\frac{1}{y-x}\right)\theta (y-x)\right]}_{+},$$

$${\int }_0^{1}dy{V}^{(0)}(x,y)G_n(y)=-\frac{{\gamma }_n^{(0)}}{2}{G}_n(x),$$

$$G_n(x)=x(1-x)C_n^{(3/2)}(2x-1).$$

Full text

1 Introduction

Rigorous quantitative understanding of the heavy-quarkonium production at high-energy colliders Bodwin:2013nua 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 Andronic:2015wma and heavy-quark Yukawa couplings to Higgs Sirunyan:2018fmm ; Aaboud:2018txb . An effective field-theoretic framework called the nonrelativistic QCD (NRQCD) Bodwin:1994jh 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$ and $\bar{Q}$ in the bound state. NRQCD is blind to the short-distance dynamics at higher energy scales of order $\gtrsim 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 ILCBehnke:2013xla , CEPCCEPCStudyGroup:2018ghi , and FCC-eeGomez-Ceballos:2013zzn 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 $\sqrt{s}$, the cross section of a quarkonium has an uncomfortably strong dependence on the large logarithm of the ratio

[TABLE]

A straightforward extrapolation of the prediction for a lower-energy process of $\sqrt{s}\lesssim 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:

[TABLE]

The resummation of such a logarithm can be made by employing the light-cone (LC) approach Lepage:1980fj ; Chernyak:1983ej or, equivalently, the soft-collinear effective theory (SCET) Bauer:2002nz . In SCET, the scattering amplitude or current-current correlator is factorized into the following factors: the hard-scattering kernel involving scales of $\sqrt{s}$, the light-cone distribution amplitude (LCDA) that represents the collinear part, and the decay constant that involves the interactions of scales $\lesssim 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, $\rho$, or $\eta$ 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 Jia:2008ep ; Jia:2010fw ; Wang:2013ywc ; Bell:2008er . 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 $\sim 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 $\eta_{c}+\gamma$ and $\eta_{b}+\gamma$. In a fixed-order perturbation theory this process was first computed in Chung:2008km at leading order (LO) and its next-to-leading order (NLO) correction was computed analytically in Sang:2009jc ; Braguta:2010mf and numerically in Li:2009ki . Up to date the $\alpha_{s}^{2}$ correction is available Feng:2015uha ; Chen:2017pyi . The relativistic correction of the order $v^{2}$ was first considered in Fan:2012dy and $\alpha_{s}v^{2}$ correction was also obtained in Xu:2014zra . The virtual $Z$ contribution was computed up to $\alpha_{s}$ correction in Chen:2013mjb ; Chen:2013itc .

The leading-logarithmic (LL) accuracy resumming $\alpha_{s}^{n}\log^{n}r$ terms was first achieved in Jia:2008ep . The quarkonium LCDAs at NLO were obtained by matching QCD onto NRQCD in Bell:2008er ; Wang:2013ywc . In Higgs or $Z$ boson decay into $J/\psi+\gamma$ processes Bodwin:2014bpa ; Bodwin:2016edd ; Bodwin:2017wdu ; Bodwin:2017pzj , the next-to-leading-logarithmic (NLL) accuracy resumming $\alpha_{s}^{n+1}\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 $\eta_{c,b}+\gamma$ 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 Jia:2018xsy . We finally summarize in Sec. 5.

2 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 Kang:2013nha ; Kang:2014qba .

[TABLE]

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 $\mu\sim\sqrt{s}$, or $\mu_{0}\sim m_{Q}$.111There is the non-perturbative scale of the order $m_{Q}v^{2}$, which is not explicitly denoted because LDMEs is not evolved in practice but rather determined at the scale $m_{Q}$ or $2m_{Q}$ as in conventional NRQCD approach. 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: $\sigma^{\text{fixed-sing}}(r;\mu)=\sigma^{\text{sing}}(r;\mu,\mu).$ The singular part of fixed-order cross section is given by

[TABLE]

where the coefficients are 222 If one expresses the $m_{P}$ in $\sigma_{0}$ in terms of $m_{Q}$ and $v^{2}$, then one finds that

$2m_{P}\left[1+\langle v^{2}\rangle c_{v^{2}}^{\textrm{sing}}\right]=4m_{Q}\sqrt{1+\langle v^{2}\rangle}\left[1+\langle v^{2}\rangle c_{v^{2}}^{\textrm{sing}}\right]\approx\left[2\sqrt{m_{Q}}\left(1-\frac{5}{12}\langle v^{2}\rangle\right)\right]^{2},$

which agrees with $d_{s}^{(0)}$ and $d_{s}^{(v^{2})}$ in Eq. (2.26) of Ref. Xu:2014zra .

[TABLE]

The born cross section $\sigma_{0}$ is given by

[TABLE]

where $\alpha$ 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, $\langle\mathcal{O}_{1}\rangle_{P}=\langle 0|\chi^{\dagger}\psi|P\rangle\langle P|\psi^{\dagger}\chi|0\rangle$ is the non-relativistically normalized LO NRQCD LDME of the production of $S$-wave pseudoscalar quarkonium $P$, and the relativistic correction to the LO NRQCD LDME, $\langle v^{2}\rangle$, is defined by

[TABLE]

where $\bm{D}=\bm{\nabla}-ig_{s}\bm{A}$ is the spartial part of the gauge-covariant derivative and $\psi^{\dagger}\overleftrightarrow{\bm{D}}\chi\equiv\psi^{\dagger}(\bm{D}\chi)-(\bm{D}\psi)^{\dagger}\chi$. $\tilde{e}_{Q}^{2}$ is given in Eq. (20) which includes the effect of both of the virtual photon and $Z$ boson propagators. To make our paper self-contained we also copy the fixed-order cross section in Sang:2009jc into App. A.

The nonsingular part is defined by subtracting the fixed-order singular part from the fixed-order cross section as

[TABLE]

where the coefficients are $c^{\text{ns}}=(c^{\text{fixed}}-c^{\text{sing}})/(-r)\approx-14.9-4.8\log r+\log^{2}r+O(r)$ and $c^{\text{ns}}_{v^{2}}=(c_{v^{2}}^{\text{fixed}}-c_{v^{2}}^{\text{sing}})/(-r)=-1/3$. Note that we pull out a prefactor $-r$ in Eq. (7) 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_{Q}^{2}/m_{Z}^{2}$ 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: $\sigma^{\text{fixed}}(r;\mu)=\sigma(r;\mu,\mu,\mu)$. Therefore, at the energies where $\log r\sim O(1)$ hence $\mu\sim\mu_{0}$, the full cross section is consistent with the fixed-order results and at higher energies where $|\log r|\gg 1$ and $\mu\gg\mu_{0}$, the resummation implemented in full cross section becomes effective. Therefore the formula in Eq. (2) 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 $\gamma^{*}$ 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

[TABLE]

where the index $I=\gamma^{*},Z^{*}$ represents the virtual bosons. The current has the vector and axial vector components: $J^{\mu}_{V}=\bar{\psi}\gamma^{\mu}\psi$ and $J^{\mu}_{A}=\bar{\psi}\gamma^{\mu}\gamma_{5}\psi$. But the axial current does not contribute due to the opposite charge conjugation and we only have $J_{V}$ in Eq. (8). The part $L_{I}^{\mu}$ contains a matrix element of initial lepton pair, a virtual boson propagator, and electroweak charges of the quarks. Its expression is

[TABLE]

where $e$ is the electromagnetic coupling, $e_{Q}$ is the fractional electric charge of the heavy quark $Q$, $\sqrt{s}$ is the CM collision energy, $\theta_{W}$ is the Weinberg angle, $g_{V}^{e}=-\frac{1}{2}+2\sin^{2}\theta_{W}$ and $g_{A}^{e}=-\frac{1}{2}$ are the vector and axial charges of electron, and $g_{V}^{Q}=T_{Q}-2e_{Q}\sin^{2}\theta_{W}$ with $T_{Q}=\pm 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. (8) can be expressed in the form factor style as

[TABLE]

where $T$ is the time-ordered product and $\epsilon_{\perp}^{\mu\nu}=\epsilon^{\mu\nu\rho\sigma}p_{\rho}k_{\sigma}/(p\cdot k)$ is the asymmetric tensor.

2.1 Light-cone distribution amplitude

The factor $G_{P}(\mu)$ is the leading-twist result in LC factorization:

[TABLE]

where $T_{H}^{P}$ is the hard-scattering kernel, $\phi_{P}$ is the LCDA of a pseudoscalar quarkonium $P$, and $f_{P}$ is the decay constant. The scale $\mu$ is an arbitrary energy scale that separates the natural scales $m_{Q}$ and $\sqrt{s}$ which are the last arguments of each functions. Each function depends on the logarithm of their ratio: the natural scale to the scale $\mu$. For simplicity, we omit the last arguments to the functions from now on.

The hard-scattering kernel $T_{H}^{P}(x;\mu)$ describes a production of quark and anti-quark pair at ${}^{1}S_{0}$ state. The one-loop expression is given in Ref. Braaten:1982yp ; Wang:2013ywc by

[TABLE]

where

[TABLE]

here $\bar{x}\equiv 1-x$. Note that $\Delta=0$ for the naive dimensional regularization (NDR) scheme Chanowitz:1979zu and $\Delta=1$ for the t’Hooft-Veltman (HV) scheme tHooft:1972tcz ; Breitenlohner:1977hr for $\gamma_{5}$ regularization.333 In the NDR scheme, $\{\gamma_{5},\gamma^{\mu}\}=0$, for the index $\mu$ in $d$ dimension, while in the HV scheme, $\gamma_{5}$ defined in 4 dimension anticommutes $\{\gamma^{\mu},\gamma_{5}\}=0$ for $\mu=0,1,2,3$ but commutes $[\gamma^{\mu},\gamma_{5}]=0$ for $\mu=4,\cdots,d$. Note that the $\delta$ in Ref. Braaten:1982yp and the $\Delta$ in Ref. Wang:2013ywc are related as $\delta=1-\Delta$. 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 $\gamma^{+}\gamma_{5}$ as

[TABLE]

where the plus components are $\gamma^{+}=\gamma^{0}+\gamma^{3}$ and $p^{+}=p^{0}+p^{3}$. $x\in[0,1]$ is the collinear momentum fraction of a quark in the quarkonium and $[z,0]$ is the gauge link that is defined by

[TABLE]

where $\mathcal{P}$ stands for path ordering, $g_{s}=\sqrt{4\pi\alpha_{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 $\phi_{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. (14) at $z=0$ and it describes a $c\bar{c}$ pair transition into a physical quarkonium.444The definition of $f_{P}$ is different from that of Ref. Wang:2013ywc by a multiplicative factor $i$. 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 Wang:2013ywc and the relativistic correction was obtained in Bodwin:2014bpa ; Wang:2017bgv . We treat $v^{2}$ and $\alpha_{s}$ corrections are of the same size and expand up to the same power. Then, the LCDA expanded up to the leading corrections in $\alpha_{s}$ and $v^{2}$ is

[TABLE]

where

[TABLE]

and the $+$ and $++$ functions are defined in App. B. The leading $\alpha_{s}$ and $v^{2}$ corrections to the decay constant $f_{P}$ are given by

[TABLE]

where $\Psi_{P}(0)$ is the wavefunction at the origin and is defined by $\Psi_{P}(0)=\langle P(p)|\psi^{\dagger}\chi|0\rangle/\sqrt{2N_{c}}$ and $|\Psi_{P}(0)|^{2}=\langle\mathcal{O}_{1}\rangle_{P}/(2N_{c})$.555We take the nonrelativistic normalization for the LDME while Ref. Wang:2013ywc takes the relativistic normalization (see Eq. (4.8)): $\langle O(^{1}S_{0})\rangle=\sqrt{2m_{P}}\langle P(p)|\psi^{\dagger}\chi|0\rangle$. The relativistic correction agrees with the result in Ref. Wang:2017bgv with $x_{0}=\bar{x}_{0}=1/2$. The renormalon ambiguity coming from the pole mass can be avoided if we replace $m_{Q}$ by $\overline{\textrm{MS}}$ mass.666See, for example, Ref. Bodwin:1998mn . Thus we replace the pole mass $m_{Q}$ in Eq. (18) with the one-loop corrected $\overline{\textrm{MS}}$ mass in Ref. Tarrach:1980up : $m_{Q}=\overline{m}_{Q}\left[1+\alpha_{s}(\overline{m}_{Q})C_{F}/\pi\right]$ and truncate higher-order contributions than our working precisions. The singular part of the cross section can be written in terms of $G_{P}(\mu)$ in Eq. (11) as

[TABLE]

where for the virtual photon $\tilde{e}_{Q}^{2}$ is $e_{Q}^{2}$, and for the virtual photon and $Z$ boson it is

[TABLE]

where $r_{Z}=m_{Z}^{2}/s$. Note that Eq. (19) is rather a fixed-order singular cross section in Eq. (3) because the functions Eqs. (12), (16), and (18) 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}(\mu)$ as

[TABLE]

2.2 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, $\mu_{0}\sim m_{Q}$ for the LCDA or $\mu\sim\sqrt{s}$ for the hard-scattering kernel $T_{H}^{P}$ to a common scale $\tilde{\mu}$, which can be chosen to be an arbitrary scale between $\mu_{0}$ and $\mu$ because the $\tilde{\mu}$ dependence should be exactly cancelled when evolutions of all the functions are combined together. One of the simple and conventional choices is to set $\tilde{\mu}=\mu$ then, the LCDA is just evolved from $m_{Q}$ to $\mu$, 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 Efremov:1979qk ; Lepage:1980fj :

[TABLE]

where the ERBL kernel $V(x,y;\alpha_{s}(\mu))$ for a pseudoscalar meson was extensively studied in the pion form factor. In the case of quarkonium the product $f_{P}\phi_{P}$ is factorized into two parts LDME and short-distance coefficients as in Eqs. (18) and (2.1) and the RG equation Eq. (22) 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 $\mu$ and the latter is from $\mu_{0}$ to $\mu$. 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 $\mu_{0}$ rather $m_{Q}v^{2}$. Then, we simply run both parts from $\mu_{0}$ to $\mu$ by using Eq. (22).

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.

[TABLE]

where the one-loop coefficient is given by

[TABLE]

and the two-loop expression can be found in Refs. Dittes:1983dy ; Sarmadi:1982yg ; Mikhailov:1984ii ; Katz:1984gf . The eigenfunction of the one-loop kernel $V^{(0)}(x,y)$ is $G_{n}$ whose eigenvalue is $-\gamma_{n}^{(0)}/2$:

[TABLE]

where $G_{n}$ is the product of the Gegenbauer polynomial $C_{n}^{(3/2)}$ and its weight Jones:2007zd :

[TABLE]

$\gamma_{n}^{(0)}$ is the LO anomalous dimension and here we follow the convention of Ref. Agaev:2010aq

[TABLE]

The NLO anomalous dimension $\gamma_{n}^{(1)}$ 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,

[TABLE]

where the $k$-th Gegenbauer coefficient of LCDA at the scale $\mu_{0}$ is

[TABLE]

with $N_{n}=4(2n+3)/[(n+1)(n+2)]$. The coefficient at LO is simple

[TABLE]

where $C^{(3/2)}_{n}(0)=(-1)^{n/2}\tfrac{(n+1)!!}{n!!}$ for even $n$ and zero for odd $n$. Explicitly, the LO ERBL equation for the $n$-th moment of the LCDA is given by

[TABLE]

Following the convention of Ref. Agaev:2010aq , the solutions of the scale evolution factor $U_{nk}(\mu,\mu_{0})$ up to NLL accuracy are given by

[TABLE]

where the value of $U_{nk}$ at LL accuracy is nonzero only for $n=k$: $\left[\alpha_{s}(\mu)/\alpha_{s}(\mu_{0})\right]^{\gamma_{n}^{(0)}/(2\beta_{0})}\delta_{nk}$. At NLL it is nonzero when $(n-k)$ is zero or, even and positive integer. $\beta_{n}$ is the beta function coefficients for $(n+1)$-th order in $\alpha_{s}$ and the explicit expressions of the two-loop anomalous dimensions $\gamma_{n}^{(1)}$ and $d_{nk}$ are copied in App. C.

Then, the RG evolved function $G_{P}(\mu)$ is given by

[TABLE]

where $\mathcal{M}^{(i,j)}$ are defined in terms of the RG evolved LCDA in Eq. (28)

[TABLE]

and

[TABLE]

is the coefficient in the expansion with Gegenbauer polynomials and it is non-vanishing only for even $n$ because $T_{H}^{P}$ 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

[TABLE]

We would like to note that the decay constant $f_{P}$ in Eq. (33) is not evolved at NLL because its one- and two-loop anomalous dimensions in Eqs. (27) and (65) are zeros: $\gamma_{0}^{(0)}=\gamma_{0}^{(1)}=0$. One can see this explicitly by taking the 0-th Gegenbauer coefficient of $\phi_{P}$ in Eq. (22), or in Eq. (31) because $\phi_{0}(\mu)=1$ for all $\mu$. We also emphasize that the relativistic correction $\mathcal{M}^{(0,v^{2})}(\mu)$ correctly resums logarithms proportional to $\langle v^{2}\rangle\alpha_{s}^{n}\log^{n}r$ by using the same RG evolution $U_{nk}$ and its expansion in $\alpha_{s}$ is given by

[TABLE]

This agrees with the logarithmic term in $\alpha_{s}\langle v^{2}\rangle$ correction Eq. (2.26) in Xu:2014zra .

Eventually we insert Eq. (33) into Eq. (19) and obtain the singular part and full cross section in Eq. (2). 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. (33), the first term $\mathcal{M}^{(0,0)}$ is computed using NLL expression of $U_{nk}$ while the other $\mathcal{M}^{(i,j)}$ terms are computed using LL expression. In the absolute square $|G_{P}(\mu)|^{2}$, we also drop higher-order terms proportional to $\alpha_{s}^{2}$ or $\alpha_{s}\langle v^{2}\rangle$, which are obtained in the product of $\mathcal{M}^{(i,j)}$ in Eq. (33) and $f_{P}$ in Eq. (18).

We also adopt the Abel-Padé method developed and used in Bodwin:2016edd ; Bodwin:2017pzj ; Bodwin:2017wdu to achieve faster numerical convergence at NLL accuracy and to deal with divergences associated with the relativistic corrections in LCDA.

2.3 $\gamma_{5}$-scheme dependence

As we can see from Eqs. (13b), (2.1), and (18), there are $\gamma_{5}$ scheme dependences in the hard kernel $T_{H}(x,\mu)$, LCDA $\phi_{P}(x,\mu)$ and decay constant $f_{P}(\mu)$, which are represented by the terms proportional $\Delta=0,1$ for NDR and HV schemes, respectively. It is easy to check that the $\Delta$ dependences of the factor $G_{P}(\mu)$ vanish at NLO without resummation or, RG evolution. However, it is not obvious whether the $\Delta$ dependences vanish or not at NLL accuracy due to additional scheme dependence that may enter in two-loop anomalous dimension $\gamma_{n}^{(1)}$. Note that $\gamma_{n}^{(0)}$ is $\Delta$-independent and so is the LL resummation. Ref. Melic:2001wb 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. Melic:2001wb , which gives

[TABLE]

Again, polynomials $G_{n}(x)$ is the eigenfunction of $\Delta V^{(1)}$ with the eigenvalue anomalous dimension

[TABLE]

where analytic expression of the anomalous dimension is given by

[TABLE]

At NLL, there are two types of $\Delta$ dependences in the amplitude $G_{P}(\mu)$. The one from NLL evolution factor $U_{nk}$ in Eq. (32) is proportional to the anomalous dimension in Eq. (40):

[TABLE]

Here $G_{P,n}^{\text{LL}}(\mu)=\left[\tfrac{\alpha_{s}(\mu)}{\alpha_{s}(\mu_{0})}\right]^{\tfrac{\gamma_{n}^{(0)}}{2\beta_{0}}}\,f_{P}^{(0)}T_{n}^{(0)}\phi^{(0)}_{n}$ is the LL amplitude, where $f_{P}^{(0)}=\sqrt{N_{c}m_{P}}\Psi_{P}(0)/m_{Q}$ is LO decay constant and $T_{n}^{(0)}$, $\phi_{n}^{(0)}$ are given in Eqs. (30) and (36).

The other type of $\Delta$ dependence is those from one-loop corrections $f_{P}^{(1)},T_{n}^{(1)},\phi_{n}^{(1)}$. The terms proportional to $\Delta$ are given by non-logarithmic constant parts

[TABLE]

Collecting three contributions above, we have

[TABLE]

In the second equality we rearranged terms into two parts, the one proportional to NLO correction evaluated at the $\alpha_{s}$ scale $\mu$ and the other proportional to the difference $\alpha_{s}(\mu_{0})-\alpha_{s}(\mu)$ which is a part of NLL resummation. The scheme independence at fixed-order NLO implies the cancellation of first part

[TABLE]

This is also confirmed by explicit computing $\Delta T_{n}^{(1)}$ and $\Delta\phi^{(1)}_{n}$ which are zero for odd $n$ and

[TABLE]

for even $n$. We would like to note an interesting relation between $\Delta T_{n}^{(1)}$ and $\Delta\gamma_{n}^{(1)}$:

[TABLE]

or, $\Delta\phi^{(1)}_{n}/\phi^{(0)}_{n}+\Delta f_{P}^{(1)}/f_{P}^{(0)}=\Delta\gamma_{n}^{(1)}/(2\beta_{0})$ equivalently. This implies that the constant term in the one-loop functions completely determines $\Delta$ dependence of two-loop evolution kernel and this ensures the cancellation of $\Delta^{(a)}_{n}$ and $\Delta^{(b)}_{n}$:

[TABLE]

where in $\Delta^{(b)}_{n}$ we eliminated $\Delta f_{n}^{(1)}$ and $\Delta\phi_{n}^{(1)}$ in favor of $\Delta T_{n}^{(1)}$ by using Eq. (46). Therefore, $\gamma_{5}$ scheme independence is valid at NLL accuracy.

At higher-order, we expect a similar pattern of cancellation between $\Delta$ dependent terms: cancellation between fixed-order terms at the same $\alpha_{s}$ scale as in Eq. (46) and cancellations between $n$-loop anomalous dimension from evolution factor and the constant terms of $(n-1)$-loop function as in Eq. (50).

2.4 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. (11). Its logarithmic structure can be schematically expressed as

[TABLE]

where $L\equiv\log(4m_{Q}^{2}/s)$. The largest logarithmic term at each $\alpha_{s}$ order is $\alpha_{s}^{n}L^{n}$ and then the next largest is $\alpha_{s}^{n}L^{n-1}$. The functions $T_{H}$ and $\phi$ are first expanded with the Gegenbauer polynomials, then each coefficient is resummed as $G_{P}^{nk}=f_{P}\,T_{n}U_{nk}\phi_{k}$, and summation over all $n$ and $k$ gives the resummed $G_{P}$. The individual $G_{P}^{nk}$ takes the following form

[TABLE]

where $C(\alpha_{s})$ is the fixed-order expansion in $\alpha_{s}$ and it does not depend on the logarithms,

[TABLE]

The fixed-order coefficient $C_{i}$ is given by fixed-order function $T_{n}$, $\phi_{k}$, $f_{P}$ and the coefficients $C_{ij}$ associated with anomalous dimensions are given by $U_{nk}$ in Eq. (32). 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. (51) are summed row-by-row, i.e., order-by-order in $\alpha_{s}$. On the other hand, in resummed perturbation theory, the series in the exponent of Eq. (52) are summed column-by-column, based on large-logarithmic power counting $L\sim 1/\alpha_{s}$. In Eq. (52), the first column is of the order $\alpha_{s}^{n}L^{n}\sim 1$ called the LL, the second column is $\alpha_{s}^{n}L^{n-1}\sim\alpha_{s}$ called the NLL, and so on. It is clear that which fixed-order terms in Eq. (53) 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. (32) is different from that of Eq. (52). For example, the non-exponent term contains logarithms: $\alpha_{s}(\mu)-\alpha_{s}(\mu_{0})=-\tfrac{\alpha_{s}(\mu_{0})^{2}}{2\pi}\beta_{0}\log(\mu/\mu_{0})+\cdots$. This is because, in Eq. (52) the second column is of $O(\alpha_{s})$ hence those terms beyond LL can be expanded and moved down from the exponent:

[TABLE]

where $\tilde{C}(\alpha_{s})$ includes two prefactors

[TABLE]

We have LL accuracy with first term in Eq. (55) and NLL with $O(\alpha_{s})$ terms in the large-logarithmic power counting $\alpha_{s}^{n+1}L^{n}\sim\alpha_{s}$. This alternative way of arranging logarithms is equivalent to Eq. (52) up to higher-order corrections than working accuracy and is the formula we use in this paper.

3 Numerical results

In this section, we list input parameters for numerical calculations then, present our results for the final state $\eta_{c,b}+\gamma$ 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.

3.1 Input parameters and NRQCD matrix elements

We use PDG values for $\overline{\text{MS}}$ mass $\overline{m}_{c}=1.275^{+0.025}_{-0.035}~{}\text{GeV}$ and $\overline{m}_{b}=4.18^{+0.04}_{-0.03}~{}\text{GeV}$, which gives the one-loop pole mass $m_{c}=1.483^{+0.029}_{-0.041}~{}\text{GeV}$ and $m_{b}=4.58^{+0.04}_{-0.03}~{}\text{GeV}$ and for Z-boson mass and width $m_{Z}=91.1876\pm 0.0021~{}\text{GeV}$ and $\Gamma_{Z}=2.4952\pm 0.0023~{}\text{GeV}$. We run the coupling constants for the electroweak using the code Global Analysis of Particle Properties (GAPP) Erler:1998sy ; Erler:1999ug and the coupling constants for the strong interaction using the 4-loop expression of the QCD beta function vanRitbergen:1997va . The CM energies of $B$-, $Z$- and Higgs factories are $\sqrt{s}=10.58,\,91.19,\,240$ GeV and the values of coupling at respective energies are $\alpha^{-1}(\sqrt{s})=130.855,\,127.916,\,127.473$, $\sin\theta_{W}(\sqrt{s})=0.233543,\,0.231201,\,0.236168$, and $\alpha_{s}(\sqrt{s})=0.1768,\,0.1184,\,0.1033$.

The NRQCD matrix elements such as the wave function at the origin and relative velocity were determined in Bodwin:2007fz by using two constraints: electromagnetic decay rate $\Gamma[\eta_{c}\to\gamma\gamma]$ of order $\alpha_{s}$ and the potential model. These values need updates due to changes in input parameters: charm-quark pole mass, scale of $\alpha_{s}$ from $m_{\eta_{c}}/2$ to $2\overline{m}_{c}$, experimental value of the decay rate. In the determinations of $\langle\mathcal{O}_{1}\rangle_{\eta_{c}}$ and $\langle v^{2}\rangle_{\eta_{c}}$, we used the same string tension $\sigma=0.1682\pm 0.0053~{}\text{GeV}^{2}$ Bodwin:2007fz and updated values for the 1-loop pole mass $m_{c}=1.483^{+0.029}_{-0.041}~{}\text{GeV}$, the mass difference between $J/\psi$ and $\psi(2S)$ $m_{2S}-m_{1S}=589.188\pm 0.028~{}\text{MeV}$, and the decay rate $\Gamma[\eta_{c}\to\gamma\gamma]=5.0\pm 0.4~{}\text{keV}$.777This is $30\%$ smaller than the value $7.2\pm 0.7\pm 2.0$ keV used in Bodwin:2007fz . Differently from Bodwin:2007fz , we do not take average with $\langle\mathcal{O}_{1}\rangle_{J/\psi}$ in the determination of $\langle\mathcal{O}_{1}\rangle_{\eta_{c}}$. In the decay rate formula, we have set the scale $\mu=2\overline{m}_{c}$. The updated values are as follow:

[TABLE]

The uncertainty includes variations of $\sigma$, $m_{c}$, $m_{2S}-m_{1S}$, $\Gamma[\eta_{c}\to\gamma\gamma]$. And we assumed the size of the neglected higher-order corrections in $\alpha_{s}$ and $v^{2}$ to be $30\%$ times the central values of $\alpha_{s}$ and $v^{2}$, respectively. Major sources of uncertainty are the variation of $\sigma$ and the assumed higher-order corrections.

We like to pay a bit more attention to using those values in Eq. (56). 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 $\eta_{c}$ we find the effect is as large as $40\%$ due to large coupling constant $\alpha_{s}(2\overline{m}_{c})\sim 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.

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 $\alpha_{s}$. This can be done by replacing the NRQCD matrix element with the experimental value of decay rate multiplied by short distance coefficient:

[TABLE]

where the experimental value for $\eta_{c}$ is $\Gamma_{\text{exp}}=5.0\pm 0.4$ keV Tanabashi:2018oca . We insert Eq. (58) into Eq. (18) and truncate higher-order terms than the working order. One of advantages using Eq. (58) is that the error propagations associated with the pole mass and the NRQCD matrix element $\langle\mathcal{O}_{1}\rangle_{\eta_{c}}$ become simpler. The pole mass $m_{Q}^{2}$ are cancelled by that of Eq. (5) in the cross section. The perturbative uncertainty obtained from scale variation of Eq. (58) largely contributes to uncertainty of $\langle\mathcal{O}_{1}\rangle_{\eta_{c}}$ and this contribution is now naturally combined in a correlated way with scale variations of the other part in the cross section. Another advantage from an empirical observation is that the 1-loop correction in the decay constant reduces significantly due to the large cancellation between $O(\alpha_{s})$ terms of the decay constant Eq. (18) and the decay rate in Eq. (58) as

[TABLE]

Note that the coefficient of $\alpha_{s}C_{F}/\pi$ reduces from $-3$ to $(8-\pi^{2})/4\approx-0.5$ and the coefficient of $\langle v^{2}\rangle$ changes from $-2$ to $-2/3$. In this way, we have better perturbative convergence between LL and NLL (LO and NLO).

3.2 Final results

Our numerical results for $\eta_{c}+\gamma$ 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

[TABLE]

The values of scales we choose for LL and NLL+NLO are $\mu_{0}=2\overline{m}_{c}$ and $\mu=\mu_{\text{ns}}=\sqrt{s}$ and for NLO we set all scales to be the same $\mu_{0}=\mu=\mu_{\text{ns}}=\sqrt{s}$. The perturbative uncertainties are estimated by varying $\mu$, $\mu_{\text{ns}}$ from its central value by a factor 2 up and down and by varying $\mu_{0}$ by a factor of $\sqrt{2}$. The uncertainties are summed in quadratures as in Kang:2013nha : $\sqrt{\delta\sigma_{\mu_{0}}^{2}+\delta\sigma_{\mu}^{2}+\delta\sigma_{\mu_{\text{ns}}}^{2}}$, where $\delta\sigma_{\mu_{i}}$ is the change of cross section by a scale variation of $\mu_{i}$. Here we do not include other sources of uncertainty to show the perturbative convergence and they will be included later in the final results in Table 1. The perturbative uncertainty (width of the band) decreases by a factor of $\alpha_{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 $\mu_{0}$: $\mu_{0}=\overline{m}_{c}$ instead of $2\,\overline{m}_{c}$. The NLL perturbative uncertainty at $\mu_{0}=\overline{m}_{c}$ tends to be asymmetric and smaller than that for $\mu_{0}=2\overline{m}_{c}$ because the lower scale variation from $\overline{m}_{c}$ by a factor of $\sqrt{2}$ moves $\mu_{0}$ close to the Landau pole and the scale dependence near this region is not monotonic. In comparison to Fig. 1 at $\mu_{0}=2\overline{m}_{c}$ we observe relatively better perturbative convergence between LL and NLL although the other is still reasonable. For these reasons we take $\mu_{0}=2\overline{m}_{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 $\eta_{b}$ is not available we use following value $\Gamma[\eta_{b}\to\gamma\gamma]=0.512^{+0.096}_{-0.094}~{}\text{keV}$ and for relative velocity $\langle v^{2}\rangle_{\eta_{b}}=-0.009^{+0.003}_{-0.003}$ taken from Chung:2010vz 888Note that there is no restriction to the positive definiteness of the matrix element $\langle v^{2}\rangle$. The matrix element $\langle v^{2}\rangle$ intrinsically contains a linear ultraviolet (UV) divergence that must be regulated. For example, if we employ dimensional regularization that is consistent with existing calulations of quarkonum decay and production rates at relative order $\alpha_{s}$ and $\alpha^{2}_{s}$, then the scaless power divergent integrals are discarded. Such subtractions of divergent contributions can lead to both positive and negative values. See Ref. Bodwin:2006dn for more details.. The central values of scales are $\mu_{0}=2\overline{m}_{b}$ and $\mu=\mu_{\text{ns}}=\sqrt{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: $\sqrt{s}=$10.58, 91.2, and 240 GeV and $Z$-boson decay branching fractions. The uncertainties in the table for charmonium ($\eta_{c}$) channel includes uncertainties of input decay rate $\Gamma^{\text{exp}}$ ($\pm 8\%$), relative velocity $\langle v^{2}\rangle$ ($\pm 30\%$) as well as perturbative uncertainties ($\pm 3\%$ or less) shown Fig. 1 and they are added in quadratures. For the bottomonium ($\eta_{b}$) the uncertainties are input decay rate ($\pm 20\%$), relative velocity ($\pm 30\%)$, perturbation ($\pm 3\%$ or less). The final uncertainties quoted in Table 1 are dominated by uncertainty of input decay rate.999We do not include relatively small uncertainties from $\overline{\text{MS}}$ mass ($\pm 2\%$) and from higher-order electroweak corrections.

There is an independent prediction for $\langle v^{2}\rangle_{\eta_{b}}=0.042$ in Ref. Li:2012rn . The authors of that reference have determined that value by making use of the Gremm-Kapustin relation. If we use this numerical value to compute the branching fraction for $e^{+}e^{-}\to\eta_{b}+\gamma$ at $\sqrt{s}=m_{Z}$ given in Table 1, we obtain $2.78\times 10^{-8}$, which is well within the prediction $(2.80\pm 0.53)\times 10^{-8}$ given in Table 1. However, we have not included the corresponding analysis into our final results listed in Table 1 because the determination of $\langle v^{2}\rangle$ by making use of either the lattice or Gremm-Kapustin approaches suffers too large uncertainties to determine even the signs of the matrix elements as is stated in Ref. Bodwin:2006dn .

4 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) Jia:2018xsy (gray line) for $\sigma(e^{+}e^{-}\to\eta_{c}+\gamma)$ at $\sqrt{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.

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 Li:2009ki and invisibly small error of NNLOChen:2017pyi only include a charm-quark 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}$Sang:2009jc , and NLO$+v^{2}$ resummationFan:2012dy use the LDME of Bodwin:2007fz , which should be updated with improved measurement of $\Gamma[\eta_{c}\to\gamma\gamma]$ as discussed around Eq. (56) and with the updated LDME, we expect decrease of the cross section by about $10\sim 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 $\mu=\sqrt{s}$ for the hard-scattering kernel and $\mu_{0}=2\overline{m}_{c}$ for the LCDA and decay constant and they are varied by a factor of 2 for $\mu$ and $\sqrt{2}$ for $\mu_{0}$ as discussed in previous section. While we use $\overline{\text{MS}}$ mass $\overline{m}_{c}=1.275^{+0.25}_{-0.35}$, many of previous results use the pole mass $m_{c}=1.4\!\sim\!1.6$ GeV. NLO Li:2009ki sets $\mu=2m_{c}$ and NLO$+v^{2}$ resummation Fan:2012dy sets $\mu=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}$ Sang:2009jc makes most conservative choice $\mu=\sqrt{s},2m_{c},$ and $m_{c}$, which leads to relatively larger uncertainty compared that of NLO$+v^{2}$ resummation. NNLO Chen:2017pyi chooses different values for renormalization scale $\mu_{r}=\sqrt{s}/2$ and the factorization scale $\mu_{\Lambda}=1.0$  GeV. We do not include the result of Braguta:2010mf because the 1-loop coefficient is not consistent with other results Sang:2009jc ; Wang:2013ywc .

Recently the Belle experiment analyzed S-wave ($\eta_{c}+\gamma$) and P-wave ($\chi_{cJ}+\gamma$ with $J=0,1,2$) channels Jia:2018xsy . While P-wave cross section ($J=1$) and upper limits ($J=0,2$) are consistent with the theoretical predictions Chung:2008km ; Sang:2009jc ; Li:2009ki ; Braguta:2010mf , S-wave upper limit $\sigma^{\text{exp}}_{\eta_{c}+\gamma}\leq 21.1$ fb at 90 % credibility level is in tension with our NLL+NLO prediction $32.7\pm{2.8}$ fb by $4.1\sigma$. This reminds us the puzzle in exclusive $J/\psi+\eta_{c}$ production Abe:2004ww ; Aubert:2005tj ; Braaten:2002fi ; Liu:2002wq , 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 Zhang:2005cha ; Bodwin:2007fz ; Bodwin:2007ga ; He:2007te .101010Recently, Feng:2019zmt reports the $K$-factor (NNLO/LO) between +20$\%$ and -40$\%$ depending on scale choice. 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 Sirunyan:2018fmm ; Aaboud:2018txb ; Bodwin:2013gca ; Mao:2019hgg . 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.

5 Summary

We resum large logarithms of $4m_{Q}^{2}/s$ at NLL accuracy for exclusive production for $\eta_{b,c}+\gamma$ 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 $\alpha_{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 $\alpha_{s}$ accuracy valid in both resummation region ($r\ll 1$) and fixed-order region ($r\sim O(1)$) where $r=4m_{Q}^{2}/s$ such that the results with the same formalism in Eq. (2) 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 $\eta_{c,b}$ is the pseudoscalar, which involves an ambiguity in handling $\gamma_{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 $\gamma_{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. (50).

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 $\alpha_{s}$ contribution that can be contained in NRQCD matrix element. By doing this all the input formula are computed at the same $\alpha_{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 $\Gamma[\eta_{b,c}\to\gamma\gamma]$ 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 differences in input parameters and uncertainty estimates. We also find Belle’s recent upper limit $21$ fb is about $4\,\sigma$ away from our prediction $33\pm 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.

Acknowledgements.

We thank Geoffrey T. Bodwin for sharing his knowledge on Abel-Padé method which was essential for resumming the non-convergent series appearing in the evolution of the heavy quarkonium LCDAs. We also thank Chaehyun Yu and Yu Jia for useful conversations during the completion of this work. The work of H.S.C. is supported by the Alexander von Humboldt Foundation. The work of D.K. is supported by NSFC through Grant No. 11875112 and the work of X.W. is supported by Fudan Scholar program. The work of J.-H.E. ​, U-R.K. ​, and J.L. is supported by the National Research Foundation of Korea (NRF) under Contract No. NRF-2017R1E1A1A01074699 (J.-H.E. ​, U-R.K. ​, J.L.), NRF-2018R1A2A3075605 (J.-H.E. ​, U-R.K.), NRF-2018R1D1A1B07047812 (U-R.K.), NRF-2019R1A6A3A01096460 (U-R.K.), and NRF-2017R1A2B4011946 (J.-H.E.). D.K. would like to thank the hospitality of QCD Group, Korea University where an important part of this work was carried out.

Appendix A fixed-order cross section

The fixed-order cross section up to $O(\alpha_{s},v^{2})$ was computed in Sang:2009jc

[TABLE]

where the coefficients $c^{\text{fixed}}$ and $c^{\text{fixed}}_{v^{2}}$ are

[TABLE]

Note that our coefficient $c^{\text{fixed}}$ is related to that of $C(r)$ in Sang:2009jc as: $c^{\text{fixed}}=3(1-r)C(r)$.

Appendix B plus distributions

Here, we give the definition of plus distributions used in the paper. The $+$ and $++$ functions are defined by

[TABLE]

The plus distribution depending on two arguments $x$ and $y$ is defined by

[TABLE]

Appendix C anomalous dimension

The NLO anomalous dimension $\gamma_{n-1}^{\parallel(1)}$ is given in Ref. GonzalezArroyo:1979df as

[TABLE]

where

[TABLE]

The off-diagonal evolution factor $d_{nk}(\mu,\mu_{0})$ is given by

[TABLE]

where

[TABLE]

and $\psi(n)$ is the digamma function.