Inclusive production of heavy quarkonium $\eta_Q$ via $Z$ boson decays within the framework of nonrelativistic QCD

In the paper, the inclusive production of heavy quarkonium $\eta_Q$ ($Q=b$ or $c$) via $Z$ boson decays within the framework of nonrelativistic QCD effective theory are studied. The contributions from the leading color-singlet and color-octet Fock states are considered. Total and differential decay widths for the inclusive decays $Z \to \eta_Q+X$ are presented. It is found that the decays $Z\to \eta_Q +X$ are dominated by the $^3S_1^{[8]}$ component, so the decays can be inversely adopted to determine the values of the long-distance matrix elements $\langle {\cal O}^{\eta_{c}}(^{3}S_{1}^{[8]})\rangle$ and $\langle {\cal O}^{\eta_{b}}(^{3}S_{1}^{[8]})\rangle$, respectively. Our numerical results show that at an $e^+e^-$ collider running at the $Z$ pole with a high luminosity around $10^{35}{\rm cm}^{-2}{\rm s}^{-1}$ (a super $Z$ factory), there are about $4.5\times 10^7$ $\eta_c$ meson events and $6.1\times 10^5$ $\eta_b$ meson events to be produced per operation year, and the inclusive decays may be used for clarifying some problems on the heavy quarkonium $\eta_Q$ and nonrelativistic QCD.


I. INTRODUCTION
Heavy quarkonia have attracted a lot of interest since the discovery of the J/ψ meson. An important reason is that they provide an ideal platform for studying the interplay between the perturbative and the nonperturbative effects in QCD. The nonrelativistic QCD (NRQCD) factorization formalism [1] provides a systematic framework to separate the short-distance and the long-distance effects in the heavy quarkonium production and decay processes. Under the NRQCD factorization, the heavy quarkonium production cross sections are expressed as the products of the short-distance coefficients (SDCs) and the long-distance matrix elements (LDMEs). The SDCs describe the production of heavy quark-antiquark pairs with proper quantum numbers, which can be calculated perturbatively. The LDMEs describe the hadronization of a produced heavy quark pair into quarkonium, which are nonperturbative in nature but can be extracted from a global fit of experimental measurements or estimated by using the QCD inspired potential models etc.
Up to now, the NRQCD factorization formalism has achieved great successes in explaining the data at the high-energy colliders [2,3]. However, there are still some challenges. For instance, the global fits of the J/ψ coloroctet (CO) LDMEs from various groups are not so consistent with each other, cf. Refs. [4][5][6][7]. Thus, it is interesting to study more quarkonium processes relating the NRQCD factorization formalism. * zhengxc@cqu.edu.cn † zhangzx@itp.ac.cn ‡ wuxg@cqu.edu.cn § hxud@cqu.edu.cn ¶ gywang@cqu.edu.cn Most studies of the quarkonia focus on the J/ψ and Υ mesons due to their high detection efficiency. For instance, the J/ψ events can be reconstructed via the decays J/ψ → l + l − (l = e, µ) with high efficiency, whose total branching ratio is ∼ 12% [8]. Contrary to the J/ψ meson, there are less studies of the η c meson production. Conventionally, the decay channel used to reconstruct the η c events is η c → γγ, and the branching ratio of this decay channel is ∼ 1.6×10 −4 [8]. Moreover, it is very difficult to record the two photons from the background in a hadron collision environment. Namely the experimental detection of the η c meson is poor. A novel proposal to reconstruct the η c events through the decay channel η c → pp has been suggested in Ref. [9], whose branching ratio is ∼ 1.5 × 10 −3 1 . This proposal opened a new way to study the η c meson at the high-energy colliders, and it has been adopted to observe the η c meson by the LHCb Collaboration [10,11]. Recent theoretical studies of the η c production at the LHC can be found in Refs. [12][13][14][15][16][17][18][19][20].
The η b meson has the same quantum numbers as those of the η c meson, but has different constituent quark mass. Since the heavier bottom quark mass, the η b meson is a better object for applying NRQCD. Thus it is interesting to study the η b and η c production applying the NRQCD factorization at the same time, although the observations on the η b are scarce. Up to now, the η b has been observed only through the feed-down contributions, i.e., from the decays of excited bottonium states. Therefore, the stud-1 In fact, the decay channels ηc → ΛΛ and ηc → ΣΣ, whose branching ratios are ∼ 1.07 × 10 −3 and ∼ 2.1 × 10 −3 , respectively [8], may also be used to identify the ηc meson so as to increase the detection efficiency of ηc. It is not very difficult to detect the strange baryon pairs produced from the ηc decay with vertex detectors because they carry high momentum from the ηc and make tracks.
ies of the η b production from various processes are requested.
At the LHC or an e + e − collider running around the Z pole and with an accessible high-luminosity (a super Z factory), the production of the heavy quarkonium η Q through Z boson decays can provide abundant information. The inclusive Z boson production cross section at the LHC with the collision energy 13 TeV is ∼ 56 nb [21]. With the luminosity of 10 34 cm −2 s −1 , there are ∼ 5.6 × 10 9 Z bosons to be produced per operation year at the LHC. A Chinese group has proposed to build a super Z factory [22], and its luminosity of the super Z factory could reach to 10 34−36 cm −2 s −1 , which is higher than that of the LEP-I by three to five orders. The Z boson production cross section is ∼ 30 nb, and there are about 3 × 10 9−11 Z bosons to be produced per operation year at the super Z factory. Therefore, it is interesting to study the production of heavy quarkonia through Z boson decays.
The production of heavy quarkonia (J/ψ and Υ etc) through Z decays has been extensively studied at the leading order in α s and v Q [23][24][25][26][27][28][29][30], the typical velocity of the heavy quark in quarkonia. For the J/ψ and Υ production through Z boson decays, the CO contributions have been estimated in Refs. [31][32][33][34], the next-to-leadingorder (NLO) QCD corrections have been calculated in Ref. [35], and the leading and next-to-leading logarithms of m Z /m Q have been resummed through the fragmentation approach in Ref. [36]. In the present paper, we devote ourselves to studying the inclusive production of η Q with Q = c or b through the Z boson decays.
According to NRQCD, for the η Q production, the leading color-singlet (CS) and color-octet (CO) Fock states are 1 S [1] 0 at v 3 Q order and 1 S [8] 0 , 3 S [8] 1 , and 1 P order. Although the CO contributions are suppressed by v 4 Q order compared to the CS contribution in the longdistance part, the CO contributions may be enhanced in the short-distance part. Therefore, besides the CS state 1 S [1] 0 , we also consider the CO states 1 S [8] 0 , 3 S [8] 1 , and 1 P [8] 1 . The remaining parts of the paper are organized as follows. In Sec.II, we briefly present the useful formulas for calculating the Z → η Q inclusive decays. In Sec.III, numerical results are presented. Section IV is reserved for discussion and conclusion.

II. CALCULATION TECHNOLOGY
Under the NRQCD factorization formalism, the decay width for the inclusive process Z → η Q + X can be written as 1 ) , and O ηQ ( 1 P [8] 1 ) are involved. To calculate the decay width for Z → η Q + X, we first calculate the decay widths for a free on shell (QQ) pair with the quantum numbers 2S+1 L ) .
The decay width for the (QQ) where N indicates that there are N particles in the final state, indicates the sum over the spin and color states of initial and final particles, and 1/3 comes from the polarization average of the initial Z boson. dΦ N is N -body differential phase space and M denotes the amplitude for the (QQ) ] pair. In the following, we sketch the formulas used in the calculation of these decay channels, successively.
A. Z → ηQ( 1 S [8] 0 , 3 S which are shown in Fig.1. The amplitude (M = M 1 + M 2 ) for the decay channel can be written down according to the two Feynman diagrams. For the 1 S [8] 0 ( 3 S [8] 1 ) case, we have where V Q and A Q are vector and axial electroweak couplings. More explicitly, V Q = T 3Q − 2e Q sin 2 θ W and A Q = T 3Q , where T 3Q and e Q are weak isospin and the charge of fermion Q in units of positron charge, respectively. p 11 = p 1 /2+q and p 12 = p 1 /2−q are the momenta of the Q andQ in the (QQ)[ 2S+1 L [8] J ] pair, and Π 1(3) is the spin-singlet (spin-triplet) projector, i.e., (12) and Λ a 8 = √ 2T a is the CO projector. For the 1 P [8] 1 case, we have Squaring the amplitude M and integrating the squared amplitude over the two-body phase space, we obtain the LO contribution for the decay channel Z → (QQ)[ 2S+1 L [8] J ] + g.

Virtual corrections
The NLO virtual corrections come from the interference of the one-loop diagrams and the LO diagrams. Four sample one-loop Feynman diagrams are shown in Fig.2. The fourth sample Feynman diagram is specific to the 3 S [8] 1 channel. The amplitude is too lengthy to be listed here.
There are UV divergences in the self-energy and vertex diagrams and IR divergences in the vertex and box diagrams. We adopt dimensional regularization with d = 4 − 2ǫ to regularize these divergences. Then the divergences appear as pole terms in ǫ. The γ 5 matrix should be noted in dimensional regularization, and we adopt the reading point prescription [37] to deal with it.
The UV divergences should be removed through renormalization. In the calculation, the renormalization scheme is adopted as follows: the renormalization of the quark field, the quark mass and the gluon field is carried out in the on-mass-shell (OS) scheme, while the renormalization of the strong coupling constant is carried out in the modified minimal subtraction (MS) scheme. The renormalization constants are where γ E is the Euler constant, µ R is the renormalization scale, β 0 = 11 − 2n f /3 is the one-loop coefficient of the QCD β function, and n f is the flavor number of active quarks. β ′ 0 = 11 − 2n lf /3 and n lf = 3 is the number of light-quark flavors. For SU (3) group, C F = 4/3, T F = 1/2, and C A = 3.
In the calculation, the threshold expansion method [38] is employed to extract the SDCs, i.e., we expand the relative momentum (q) of the (QQ)[ 2S+1 L [8] J ] pair before performing the loop integration. Then the Coulomb divergences, which are IR power divergences and vanish in dimensional regularization, do not appear in our calculation.

Real corrections
The real corrections to the decay channel Z → (QQ)[ 2S+1 L [8] J ] + g come from the processes Z → (QQ)[ 2S+1 L [8] J ] + gg and Z → (QQ)[ 2S+1 L [8] J ] + qq, where q = u, d, s. In the calculation, we use i ǫ µ * i ǫ ν i → −g µν to sum the polarizations of the final-state gluons. The unphysical polarization contributions are subtracted through the process involving ghost-pair production, i.e., Z → (QQ)[ 2S+1 L [8] J ] + u gūg . Six sample Feynman diagrams for the real corrections are shown in Fig.3. The sixth diagram is specific to the 3 S [8] 1 channel. There are IR divergences in the real corrections. These IR divergences should be regularized by dimensional regularization as those in the virtual corrections. In order to simplify the calculation of the real corrections under dimensional regularization, we adopt the two-cutoff phasespace slicing method [39] to isolate the divergent terms. Under this method, the differential decay width for the real corrections can be decomposed into three parts, where dΓ S denotes the contribution from the phase space In the calculation, the two cutoff parameters should be taken as δ c ≪ δ s ≪ 1. Applying the eikonal and collinear approximations to the soft and hard-collinear parts respectively, Γ S and Γ HC can be calculated analytically in d space-time dimensions. Due to the constraints Applying the NRQCD factorization to the production of an on shell (QQ)[ 1 P [8] 1 ] pair, the decay width can be written as dΓ 1 P [8] 1 =dΓ 1 P [8] 1 O (QQ)[ 1 P [8] 1 ] ( 1 P [8] 1 ) where dΓ 1 P [8] 1 denotes the decay width for an on shell (QQ) pair with quantum numbers 1 P [8] 1 , O (QQ)[ 1 P [8] 1 ] ( 2S+1 L [8] J ) denotes the LDMEs for the quark pair. The LDME O (QQ)[ 1 P [8] 1 ] ( 1 P [8] 1 ) starts at α 0 0 ) start at α s order. Up to α 2 s order, the second term vanishes in this decay channel. Then, we have .
For the decay channel Z → (QQ)[ 1 S [8] 0 ] + Q ′Q′ , there are two Feynman diagrams which are shown in Fig.7. The amplitude can be written as M = M 1 + M 2 . The amplitudes M 1 and M 2 are the same as M 5 and M 6 in Eqs. (23) and (24), but here we have  (24) through replacement Π 1 → Π 3 . The amplitudes M 3 and M 4 , which correspond to the two diagrams in Fig.8, are as follows: For the decay channel Z → (QQ)[ 1 P [8] 1 ] + Q ′Q′ , there are two Feynman diagrams which are shown in Fig.7. The amplitude can be written as M = M 1 + M 2 . The amplitudes M 1 and M 2 have the same form as M 5 and M 6 in Eqs. (31) and (32), but here p 2 There are six Feynman diagrams for the decay channel Z → (QQ)[ 1 S [1] 0 ] + gg. Half of the Feynman diagrams are shown in Fig.9, and the other three Feynman diagrams can be obtained from these diagrams through the substitution p 2 ↔ p 3 . According to those diagrams, the amplitude (M = 6 i=1 M i ) of the process can be written down, and we have The other three amplitudes M i (i = 4, 5, 6) can be obtained from M i (i = 1, 2, 3) via the substitution p 2 ↔ p 3 .

III. NUMERICAL RESULTS
In the calculations, the package FeynArts [40] is employed to generate Feynman diagrams and amplitudes, the package FeynCalc [41,42] is employed to carry out the color and Dirac traces, the package $Apart [43] is employed to conduct the partial fraction, the package FIRE [44] is employed to do the integration-by-parts (IBP) reduction, and the package LoopTools [45] is used to compute the one-loop master integrals numerically. The phase-space integrations are performed by using the package Vegas [46].
The necessary input parameters for the numerical calculation are taken as follows: m c = 1.5 GeV, m b = 4.75 GeV, m Z = 91.1876 GeV, sin 2 θ W = 0.231, α = 1/128, (38) where α is the electromagnetic coupling constant at m Z . For the strong coupling constant, we adopt the one-loop formula .
According to α s (m Z ) = 0.1179 [47], we obtain α s (2m c ) = 0.234 and α s (2m b ) = 0.175. For the LDMEs, we derive the LDMEs for the η Q from the experimentally extracted LDMEs for the J/ψ(Υ) via the heavy quark spin symmetry (HQSS), i.e., These relations are expected to hold to relative order v 2 Q . Several sets of the LDMEs for the J/ψ and the Υ extracted from the global fits by several groups are listed in Tables I and II. The factorization scale of the LDMEs has been taken as µ Λ = 1.5 GeV for the J/ψ and the Υ in Tables I and II. Thus, we also take this value for µ Λ in this paper.

A. Integrated decay widths
In this subsection, we give the decay widths for different decay channels and the total decay widths for the inclusive η Q production via Z boson decays.
The decay widths for the decay channels contributing to Z → η c + X are given in Tables III, IV, V and  LDMEs Gong Feng et al. [48]   J ) + g based on three sets of LDMEs, where "NLO" denotes the results up to NLO accuracy. The very large NLO correction in the 3 S [8] 1 case comes from the contribution of the real correction processes Z → ηc( 3 S [8] 1 ) + qq. Table III, the decay widths for the decay channels Z → η c ( 2S+1 L [8] J ) + g up to LO and NLO accuracy in α s are presented. We can see that the NLO correction is larger than the LO contribution in the 1 S 1 cases. The reason of the large NLO correction in the 1 S [8] 0 case is that only the vector coupling of the Z − cc vertex contributes to the decay width of Z → η c ( 1 S [8] 0 ) + g at the LO level, while both the vector and axial-vector couplings of the Z − cc vertex contribute to the NLO correction through the real corrections. Moreover, the strength of the axial-vector coupling is stronger than that of the vector coupling in the Z − cc vertex. Thus, the large NLO correction in the 1 S [8] 0 case is expected. The reason for the very large NLO correction in the 3 S [8] 1 case is that there are gluon fragmentation diagrams for the processes Z → η c ( 3 S [8] 1 ) + qq at αα 2 s order. One of the gluon fragmentation diagram is the sixth diagram in Fig.3. The decay width for Z → η c ( 3 S [8] 1 ) + qq(q = u, d, s) 2 is 1.15 × 10 3 keV under the LDME extracted by Chao et al, which is very close to the decay width for Z → η c ( 3 S [8] 1 ) + g at the NLO level.

VI. In
From these tables, we can see that the dominant contributions come from the decay channels Z → η c ( 1 S These dominant decay channels can be understood by the fragmentation mechanism. For the decay channel Z → η c ( 1 S [1] 0 ) + cc, the decay width is dominated by the (anti)quark fragmentation process of Z → cc followed by c(c) → η c . In this fragmentation process, the quark propagator is of order 1/m ηc , and the gluon propagator is of order 1/m 2 ηc . For the decay channels Z → η c ( 1 S [8] 0 , 1 P [8] 1 ) + cc, the decay widths are also dominated by the (anti)quark fragmentation process. However, the involved LDMEs in the two decay channels are suppressed by powers of v c compared with the CS LDME. Thus, the decay widths of the decay channels Z → η c ( 1 S [8] 0 , 1 P [8] 1 ) + cc are suppressed compared with that of Z → η c ( 1 S [1] 0 ) + cc. For the decay channels 2 When we present the results for Z → η Q ( 3 S [8] 1 ) + qq individually, we actually give the contribution from the fragmentation diagrams, which is gauge invariant and counts almost the whole contribution of the NLO decay width of Z → ηc( 3 S [8] 1 ) + g. Z → η c ( 3 S [8] 1 ) + cc(bb, qq), the decay widths are dominated by the fragmentation processes of Z → cc(bb, qq) followed by a quark or an antiquark fragments into the η c and Z → ccg(bbg, qqg) followed by g → η c . In the fragmentation processes Z → cc(bb, qq) followed by c(b, q,c,b,q) → η c , the quark propagator is of order 1/m ηc , and the gluon propagator is 1/m 2 ηc . Since the gluon propagator is fixed as 1/m 2 ηc in the whole phase space, the CO channels Z → η c ( 3 S [8] 1 ) + cc(bb, qq) are more important than the CS channel Z → η c ( 1 S [1] 0 ) + cc although the CO channels are suppressed by powers of v c compared with the CS channel. The decay channels Z → η c ( 1 S [8] 0 , 3 S [8] 1 , 1 P [8] 1 ) + g at LO in α s have no fragmentation contribution, the quark propagator in these channels is of order 1/m Z . Other channels Z → η c ( 1 S [8] 0 , 1 P [8] 1 )+bb and Z → η c ( 1 S [1] 0 ) + gg also have no fragmentation contribution, thus they are suppressed. 3 Due to the fact that the SDCs of the 3 S [8] 1 channels are greatly enhanced compared to other channels and the 3 S [8] 1 channels dominate the decay Z → η c + X, the total decay width of Z → η c + X is sensitive to the CO LDME O ηc ( 3 S [8] 1 ) . Therefore, the process Z → η c +X provides a good platform to determine the value of O ηc ( 3 S [8] 1 ) . Moreover, according to HQSS, we have O J/ψ ( 1 S [8]  Summing the contributions from the considered decay channels, we obtain the decay width for the inclusive process Z → η c + X which is given in Table VII.  The contributions to the decay width of Z → η b + X from the considered decay channels are given in Tables VIII, IX, X and XI. The very large NLO correction to the decay channel Z → η b ( 3 S [8] 1 )+g comes from the processes Z → η b ( 3 S [8] 1 ) + qq whose decay width is 8.21 keV under the LDME extracted by Gong 1 ) + cc and Z → η b ( 3 S [8] 1 ) + qq due to the fragmentation mechanism in these channels. Different from the η c case, among these dominant decay channels, the CS channel Z → η b ( 1 S [1] 0 ) + bb is the most important channel. There are two reasons: One is that m Z /m η b is smaller than m Z /m ηc , which leads to the enhancement in the SDCs of the 3 S [8] 1 channels is weakened for the η b case; The other is that the CO LDME O ηQ ( 3 S [8] 1 ) is more suppressed compared with the CS LDME O ηQ ( 1 S [1] 0 ) in the η b case. However, for the 3 S 1 channels dominate the decay Z → η b + X. Therefore, the decay Z → η b + X can be used to determine the value of O η b ( 3 S  Summing the contributions from the considered decay channels, we obtain the decay width for the inclusive process Z → η b + X which is presented in Table XII.

B. Differential decay widths
In this subsection, we present the differential decay widths dΓ/dz for Z → η Q + X, where the energy fraction is defined as z ≡ 2p ηQ · p Z /p 2 Z . Since the dominant contributions to the decay process Z → η Q + X come from the 1 S 1 channels in this subsection. The differential decay widths dΓ/dz for Z → η c + X based on the three sets of LDMEs are given in Fig.10. From the figure, we can see that the 3 S [8] 1 channels dominate the decay Z → η c + X, which was also shown in the last subsection by the integrated decay widths. The distributions of the CS channel and the CO channels have different shapes. The curve of the CS channel has a peak at a moderate z value, while the curves of the CO channels have a peak at a small z value. This feature can be used to determine the CS LDME O ηc ( 1 S [1] 0 ) and the CO LDME O ηc ( 3 S [8] 1 ) more precisely. The differential decay widths for Z → η b + X based on the two sets of the LDMEs are given in Fig.11. Here, the differential decay width of the channel η b ( 1 S [1] 0 ) + gg is also given. Similar to the η c case, the curves of the CS channels have a peak at a moderate z value while the curves of the CO channels have a peak at a small z value. However, since the CS contribution is comparable with the CO contribution in the η b case, the shape of the inclusive process Z → η b + X is significantly different from that of Z → η c + X.  [4], the middle one shows the distributions based the LDMEs of Chao et al [5], and the bottom one shows the distributions based the LDMEs of Gong et al [6].

IV. DISCUSSION AND CONCLUSION
In the present paper, we have studied the inclusive production of η Q (Q=c or b) through Z boson decays. The CS ( 1 S [1] 0 ) and the CO ( 1 S [8] 0 , 3 S [8] 1 , and 1 P [8] 1 ) Fock states are considered. The integrated and differential decay widths for the related channels are computed, and FIG. 11. The differential decay widths dΓ/dz for Z → ηc +X, where µR = 2m b . The top one shows the distributions based on the LDMEs of Gong et al [48], and the bottom one shows the distributions based the LDMEs of Feng et al [49].
the results show that the decay width of Z → η Q + X is dominated by the CO 3 S [8] 1 production. It means that the decay width of Z → η Q + X is sensitive to the value of the LDME O ηQ ( 3 S [8] 1 ) . Hence, the two processes Z → η c + X and Z → η b + X can be used to determine the values of O ηc ( 3 S [8] 1 ) and O η b ( 3 S [8] 1 ) . Moreover, via HQSS, the measured value of O ηQ ( 3 S [8] 1 ) from the process Z → η Q + X can also give a certain constraint on the value of O J/ψ(Υ) ( 1 S [8] 0 ) . Note that this conclusion depends partly on the exact values of LDMEs, and it is applicable only if the LDMEs are at the same order of magnitude of the ones quoted in this paper.
The differential distributions dΓ/dz are shown in fig-ures. The distributions of the CS and the CO components are very different. The distributions of the CS component have a peak at a moderate z value, while the distributions of the CO components have a peak at a small z value. Thus, the CS LDME O ηQ ( 1 S [1] 0 ) and the CO LDME O ηQ ( 3 S [8] 1 ) can be determined more precisely through measuring the energy distribution of the process Z → η Q + X.
In a hadronic collider such as the LHC, the production of the heavy quarkonium η Q is dominated by hadronic production, so it is difficult to pick up the production events via Z decays. Thus the calculations on the production via Z boson decays here may be really useful as reference mainly for the production in a super Z factory. The total cross section for the η Q production via the electron and positron annihilation at the Z pole, e + e − → Z → η Q + X 4 can be derived from the decay width Γ Z→ηQ+X through the formula derived in the Appendix A1 of Ref. [56], i.e., σ e + e − →ηQ+X = e 2 (1 − 4sin 2 θ W + 8sin 4 θ W ) Then we obtain σ e + e − →ηc+X = 45.0 pb, (41) σ e + e − →η b +X = 0.608 pb, (42) where the input values for the LDMEs have been taken as those extracted by Chao et al [5] and Gong et al [48], which have been presented in Tables I and II. If the luminosity of a Z factory can be up to 10 35 cm −2 s −1 [22], then there are about 4.5 × 10 7 η c and 6.1 × 10 5 η b to be produced per operation year. Therefore, at a high luminosity Z factory with highly rejecting backgrounds, those two production processes can be studied thoroughly.