Next-to-leading order QCD corrections to $Z\to \eta_Q+Q+\bar{Q}$

It has been found that at a high luminosity $e^+ e^-$ collider, sizable $\eta_c+c\bar{c}X$ and $\eta_b+b\bar{b}X$ events can be produced when it works around the $Z$ peak. In this paper, we calculate the decay widths of $Z \to \eta_c+c+\bar{c}+X$ and $Z \to \eta_b+b+\bar{b}+X$ up to next-to-leading order (NLO) accuracy. We find that the NLO corrections are significant in these two processes. After including the NLO corrections, the decay widths of $Z \to \eta_c+c+\bar{c}+X$ and $Z \to \eta_b+b+\bar{b}+X$ are enhanced by about $37\%$ and $28\%$ for the case of $\mu_R=2m_c$ and $\mu_R=2m_b$, respectively. The differential decay widths $d\Gamma/ds_1$ and $d\Gamma/dz$ for these two decay processes are also analyzed.


I. INTRODUCTION
Heavy quarkonium production presents an ideal laboratory for the study of the interplay between the perturbative and nonperturbative effects of QCD; it has been a focus of theoretical and experimental interest since the discovery of J/ψ in 1974.In order to describe the quarkonium production, the color-evaporation model (CEM) [1,2], the color-singlet model (CSM) [3][4][5], and the nonrelativistic QCD (NRQCD) effective theory [6] have been proposed.Among them, the NRQCD effective theory provides a systematic way of separating the short-distance and long-distance effects in the quarkonium production, and has achieved great success in describing the experimental data of the quarkonium production, especially for the unpolarized cross section of the J/ψ hadroproduction [7][8][9][10].However, there are still challenges to NRQCD.For instance, the hadroproduction cross section of η c measured by the LHCb experiments [11] can be well described by the color-singlet contribution, i.e., the color-octet contribution should be very small [12].This seems inconsistent with the heavyquark spin symmetry (HQSS) relation between the longdistance matrix elements (LDMEs) of η c and J/ψ. 1 It is important to study more processes involving the charmonium for testing the NRQCD factorization.
It has been found that the heavy quarkonium production through Z boson decays can provide a good platform for studying the production mechanism of quarkonia, which has attracted great attention .A large number of Z boson events can be accumulated at the LHC or a future high-luminosity e + e − collider running around the Z pole.At the LHC, there are about 10 9 Z bosons to be produced per year [29].It is well known that several proposed high-luminosity e + e − colliders, such as the ILC [40], FCC-ee [41], CEPC [42], and Super Z Factory [43], are planned to run at the Z pole for a period of time.When the e + e − collider runs at the Z pole and with a luminosity of 10 34−36 cm −2 s −1 , there are about 10 9−11 Z bosons to be produced per year [44].These colliders will open new opportunities for studying the quarkonium production through Z boson decays.
Most studies on the heavy quarkonium production through the Z boson decays focus on the spin-triplet J/ψ and Υ production, while few studies are for the spinsinglet η Q (Q = b or c) production.In our recent work [38], we studied the inclusive production of the η Q via the Z boson decays up to order αα 2 s within the framework of NRQCD, in which the leading color-singlet ( 1 S [1] 0 ) and color-octet 1 , and 1 P [8] 1 ) Fock states are considered.The study found that there are many interesting features in these production processes.An important channel contributing to the inclusive production Experimentally, its decay width can be measured separately through the heavy-flavor tagging technology.Therefore, it is helpful to do a precise theoretical study on this channel.In this paper, we devote ourselves to studying the decay Z → η Q + Q + Q + X, which starts at order αα 2 s , up to NLO QCD accuracy.We will use the CSM, which is the leading-order (LO) contribution (in v Q , where v Q is the velocity of the heavy quark or the heavy antiquark in the quarkonium rest frame, v 2 c ≈ 30% for the η c and v 2 b ≈ 10% for the η b [45]) of NRQCD, 2 to calculate the decay width of 2 The next-order relativistic correction to the color-singlet contribution is suppressed by order v 2 Q , while the color-octet contribution is suppressed by order v 4 Q .It is noted that the shortdistance factor of the color-octet contribution may be enhanced compared to that of the color-singlet contribution.In this work, we assume the color-octet contribution is very small, and focus on the color-singlet contribution.
The NLO QCD corrections to Z → η Q + gg have recently been finished through the CSM [39].The authors there found that the NLO corrections are significant due to the fragmentation diagrams appearing at the NLO level. 3Reference [39] and the present paper give a complete study on the η Q production through Z boson decays up to NLO QCD accuracy under the CSM.
The remaining parts of the paper are organized as follows.In Sec.II, we briefly present useful formulas for the process Z → η Q + Q + Q + X at the LO accuracy.In Sec.III, we present the formulas for calculating the NLO QCD corrections to the process Z → η Q + Q + Q + X.In Sec.IV, numerical results and discussions are presented.Section V is reserved as a summary.

II. LO DECAY WIDTH
Under the NRQCD factorization, the decay width for where d Γ are short-distance coefficients (SDCs) and O ηc (n) are LDMEs.The sum extends over all of the intermediate color-singlet and color-octet states 2S+1 L [1,8] J .In the lowest-order nonrelativistic approximation (i.e., the CSM), only the color-singlet contribution with n = 1 S [1] 0 needs to be considered.In the practical calculation, we first calculate the decay width for a free on-shell (Q Q) pair with quantum number 0 ]+Q+ Q+X .Then the decay width for the η Q meson can be ob- 0 ) . 3 The large fragmentation contribution in the NLO corrections of Z → η Q + gg can be calculated through the fragmentationfunction approach, where the large logarithms of m Z /m Q in higher-order corrections can be resummed through the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution of the fragmentation functions [46].
At the LO level, there are four Feynman diagrams for the decay process Z → (Q Q)[ 1 S [1] 0 ] + Q + Q, which are shown in Fig. 1.Corresponding to the four Feynman diagrams, the LO amplitude for this process can be written as Here, V Q and A Q are vector and axial electroweak couplings, respectively.More explicitly, (6) and Λ 1 is the color projector for color-singlet state where 1 is the unit matrix of the SU (3) c group.With these amplitudes, the LO decay width for the (Q Q)[ 1 S [1] 0 ] pair can be calculated through dΓ where denotes the sum over the spin and color states of initial and final particles.dΦ 3 is the differential phase space for the three-body final state, and where d is the number of the space-time dimensions.With these formulas, the LO decay width for Z →

III. NLO CORRECTIONS
At the NLO level, the virtual and real corrections need to be calculated.There are ultraviolet (UV) and infrared (IR) divergences in virtual correction, and IR divergence in real correction.The conventional dimensional regularization with d = 4 − 2ǫ is employed to regularize both UV and IR divergences throughout this paper.In dimensional regularization, the γ 5 problem is notorious, and we adopt a practical prescription proposed in Ref. [47].In the following subsections, we explain our main steps in calculating the virtual and real corrections.
The virtual correction at the NLO level comes from the interference of the one-loop Feynman diagrams and the LO Feynman diagrams.Four typical one-loop Feynman diagrams are shown in Fig. 2. It is noted that, compared to the J/ψ case [28], there are new type Feynman diagrams, in which the (Q Q)[ 1 S [1] 0 ] pair is produced from two virtual gluons, need to be calculated in the η Q case.These new type Feynman diagrams do not contribute to the J/ψ case.One of the new type diagrams is shown by the fourth diagram in Fig. 2.
The virtual correction to the decay width of the process where M Virtual is the amplitude for the virtual correction, and dΦ 3 is the differential three-body phase space, which has been presented in Eq.( 9).In order to take the lowest-order nonrelativistic approximation, we need to expand the amplitude in q (the relative momentum between the quark and antiquark in the (Q Q)[ 1 S [1] 0 ] pair).In the actual calculation, we expand the amplitude in q (it is equivalent to taking q = 0 here) before performing the loop integration.In the language of method of region [48], this amounts to directly calculating the contributions from the hard region.The Coulomb divergence, which is power IR divergence, will vanish in the calculation under dimensional regularization.
There are UV and IR divergences in the loop integrals.The IR divergences from the virtual correction will be canceled by the IR divergences from the real correction.The UV divergences need to be removed through renormalization.The renormalization scheme is taken as follows: For the renormalization of the heavy quark field, the heavy quark mass and the gluon field, the onmass-shell (OS) scheme is adopted, while for the renormalization of the strong coupling constant, the modified minimal subtraction (MS) scheme is adopted.With this renormalization scheme, the quantities δZ i ≡ Z i − 1 can be derived [49] where µ R is the renormalization scale, γ E is the Euler constant.
, we consider the charm-quark loop in the gluon self-energy but neglect the bottom-quark and top-quark loops in the gluon self-energy, i.e. n f = n lf +1 = 4; when µ R ∈ [m b , m t ), we consider the charm-quark and bottomquark loops in the gluon self-energy but neglect the topquark loop in the gluon self-energy, i.e. n f = n lf + 2 = 5.For SU (3) c group, In the calculations, the package FeynArts [50] is employed to generate Feynman diagrams and amplitudes, the package FeynCalc [51,52] is employed to carry out the color and Dirac traces, the packages $Apart [53] and FIRE [54] are employed to conduct partial fraction and integration-by-parts (IBP) reduction.After the IBP reduction, all one-loop integrals are reduced into master integrals, and the master integrals are calculated by the package LoopTools [55] numerically.The final phasespace integrations are calculated with the help of the package Vegas [56].
FIG. 3. Four typical real-correction Feynman diagrams for the decay process, The real correction to the process ).Four typical Feynman diagrams are shown in Fig. 3. Compared to the J/ψ case, we need to deal with new type Feynman diagrams, in which the 0 ] pair is produced from the gluon fragmentation, such as the fourth diagram of Fig. 3.
Using these Feynman diagrams, the amplitude (M Real ) for the real correction can be written down directly.Then the differential decay width for the real correction can be calculated through dΓ where dΦ 4 is the differential four-body phase space, There are IR divergences in the real correction, which come from the phase-space integration over the region where the momentum of the final gluon is close to zero.We employ the two-cutoff phase-space slicing method [57] to isolate the IR divergences in the real correction.Due to the fact that there is no collinear divergence in the present process, we only need to introduce one cutoff parameter δ s .Then the phase space for the real correction is divided into two regions: The soft region with E 4 ≤ m Z δ s /2 and the hard region with E 4 > m Z δ s /2.Here, we define the gluon energy E 4 in the rest frame of the initial Z boson.More explicitly, the real correction can be divided into two parts dΓ where and Applying the eikonal approximation to the amplitude in the soft region [57,58], we obtain Up to O(δ s ) corrections, the differential phase space for the soft region can be factorized as [57] where dΦ 3 denotes the differential three-body phase space without emitting a gluon, whose expression has been shown in Eq.( 9).Inserting Eqs.( 17) and ( 18) into Eq.(15), and carrying out the integration over p 4 [59,60], we obtain dΓ where , where E 2 and E 3 are also defined in the rest frame of the initial Z boson.Due to the constraint E 4 > m Z δ s /2 for the hard region, the contribution from the hard region is finite, then 0 ] H can be numerically calculated in four dimensions.The real correction can be obtained by summing the contributions from the hard and soft regions easily.Both the contributions from the soft and hard regions are separately dependent on the cutoff parameter δ s , while the sum of these two contributions should be independent to the choice of δ s (δ s should be small enough.).Verifying this δ s independence is an important test of the correctness of the calculation.We have checked the δ s independence, and have found that the results are independent of δ s within the error of the numerical integration when δ s varies from 10 −5 to 10 −7 .
The net NLO corrections can be obtained through summing the virtual and real corrections.After summing the virtual and real corrections, the UV and IR divergences are exactly canceled, and the finite results are obtained.The decay width dΓ Z→ηQ+Q QX can be obtained from dΓ Z→(Q Q)[ 1 S [1] 0 ]+Q QX by multiplying a factor , where R ηQ (0) is η Q radial wave function at the origin, which can be calculated by using the potential model [61].

IV. NUMERICAL RESULTS AND DISCUSSIONS
The input parameters for the numerical calculation are taken as follows [62]: m c = 1.67 ± 0.07 GeV, m b = 4.78 ± 0.06 GeV, m Z = 91.1876GeV, sin 2 θ W = 0.231, α = 1/128, (20) where m c and m b are the pole masses, α is the electromagnetic coupling constant at m Z .For the running strong coupling constant, we use the two-loop formula where is the two-loop coefficient of the QCD β function.According to α s (m Z ) = 0.118 [62], we obtain Λ With the values for Λ QCD , the strong coupling constant at any scale can be directly calculated through Eq.( 21).For the radial wave functions at the origin, we adopt the values from the potentialmodel calculation [61], i.e.,

A. Integrated decay widths
In this subsection, we give the integrated decay widths for the decay channel  In order to have a glance on the size of the NLO corrections, we first present the decay widths when the input quark masses are taken as their central values (i.e.,  In Figs. 4 and 5, the dependence of the decay widths on the renormalization scale is shown.After including the NLO corrections, the dependence of the decay widths on the renormalization scale is weakened.For Z → η c + c + c + X(Z → η b + b + b + X), the decay width decreases by 77%(57%) at LO, and by 63%(39%) at NLO when µ R varies from 2m Q to m Z .However, this dependence is still strong even including the NLO corrections.Now, let us estimate the theoretical uncertainties for these decay widths.The main uncertainty sources include the renormalization scale, the heavy quark masses and the radial wave functions at the origin4 .For the uncertainties caused by the renormalization scale, we estimate them by varying the renormalization scale between two physical energy scales involved in the processes, i.e. 2m Q and m Z .Furthermore, we take the average values of the decay widths under the two choices of the renormalization scale as their central values.For the uncertainties caused by the heavy quark masses, we estimate them by varying the heavy quark masses in the ranges given in Eq.( 20), i.e., m c = 1.67 ± 0.07 GeV and m b = 4.78 ± 0.06 GeV.For the radial wave functions at the origin, the authors of Ref. [61] did not give an error estimate.Since the potential used in Ref. [61] does not include the spin effect, the wave functions calculated in this way are accurate up to corrections of relative order v 2 Q .Therefore, we estimate the uncertainties by attaching an error of 30% of the central value for η c , and 10% of the central value for η b .More explicitly, we take |R ηc (0)| 2 = 0.810 ± 0.243 GeV Here, the first error is caused by the renormalization scale, the second error is caused by the heavy quark mass, and the last error is caused by the radial wave function at the origin.From Eqs.( 23) and ( 24), we can see that the largest error arises from the renormalization scale uncertainty for both η c and η b cases.Furthermore, we find that although the K factors are sensitive to the renormalization scale, they are insensitive to the heavy quark mass, e.g., when we vary the charm (bottom) quark mass from 1.60GeV (4.72GeV) to 1.74GeV (4.84GeV), the K factor changes from 1.50 (1.43) to 1.49 (1.44) for the η c (η b ) case.Adding the uncertainties in quadrature, we obtain  The differential distributions contain more information than the integrated decay widths, which can be used to test the current theory.Therefore, it is interesting to see the differential distributions of the two Z boson decay processes.The energy fraction carrying by η c or η b in the processes can be defined as z ≡ 2 p 0 • p 1 /m 2 Z .The LO and NLO differential decay widths dΓ/dz for the processes Z → η c +c+c+X and Z → η b +b+ b+X are shown in Figs. 6 and 7, respectively.Figs. 6 and 7 confirm the importance of the NLO corrections.For the η c production, the magnitude of dΓ/dz is increased obviously at small and moderate z values and decreased slightly at higher z values.And for the η b production, the magnitude of dΓ/dz is increased at all z values.The uncertainties for dΓ/dz are also shown in the figures, which are obtained by combining the uncertainties of the renormalization scale, the heavy quark mass and the wave function at the origin.The momenta of heavy quarks in the final state can be determined experimentally using the heavy-flavor tagging technology.Therefore, the differential decay widths dΓ/dm 12 and dΓ/dm 23 can be measured experimentally, where m 12 ≡ (p 1 + p 2 ) 2 and m 23 ≡ (p 2 + p 3 ) 2 are invariant masses of two final-state particles.We present the LO and NLO differential decay widths dΓ/dm 12 and dΓ/dm 23 for Z → η c + c + c + X and Z → η b + b + b + X in Figs. 8, 9, 10 and 11, respectively.The uncertainties for these differential decay widths are also shown in these figures.From the figures, we can see that the differential decay widths dΓ/dm 12 and dΓ/dm 23 are changed significantly after including the NLO corrections, especially for the decay Z → η c + c + c + X.Moreover, it is found that the NLO differential decay widths dΓ/dm 12 and dΓ/dm 23 are negative at the large m 12 or m 23 for Z → η c +c+c+X.This indicates that the NLO corrections are negative and larger than LO contribution in these phase space regions.In the boundary regions of the phase space, some large logarithmic terms often appear in the coefficients of the perturbation expansion, which may spoil the convergence of the perturbation series.In these boundary regions, it is necessary to resum these large logarithmic terms so as to obtain precise theoretical results.Because of these large logarithmic terms, if we calculate to a certain order (e.g.NLO) in these regions, we may obtain nonphysical negative results.Fortunately, in the considered processes, the absolute values of the negative contributions of these regions are not large.Thus, we believe that the differential decay widths for most regions of the phase space and the integrated decay widths, which are obtained in this paper, are reliable.

V. SUMMARY
In the present paper, we have studied the decays Z → η c + c + c + X and Z → η b + b + b + X up to NLO QCD accuracy.Integrated and differential decay widths of both decay processes are obtained, and the uncertainties for them are estimated.We find that the NLO corrections to the decay widths for Z → η c + c + c + X and Z → η b + b + b + X are significant.The dependence of the decay widths on the renormalization scale is very strong although the dependence is weakened after including the NLO corrections.This brings a big uncertainty to the theoretical predictions under the conventional renormalization scale setting.The higher-order corrections can reduce the uncertainty caused by the renormalization scale.However, it is very difficult to calculate the higherorder corrections beyond the NLO for these processes at present.In the literature, the principle of maximum conformality (PMC) scale-setting approach [63][64][65][66] has been suggested to eliminate such scale uncertainty, whose key idea is to determine the correct momentum flow of the process by using the nonconformal β-terms that govern the α s running behavior with the help of renormalization group equation.As a comparison, we also calculate the integrated decay widths under the PMC scale setting.Following the standard PMC procedures to the decay width of Z → η Q + Q + Q + X, we obtain Γ PMC Z→ηc+ccX = 95.4 +34.1 −32.2 keV and Γ PMC Z→η b +b bX = 14.7 +1.7 −1.7 keV.Here, the PMC predictions of the decay width are independent to the choices of µ R , and the errors come from the uncertainties of the heavy quark masses and the wave functions at the origin.
The cross sections for e + e − → Z → η Q + Q + Q + X at the Z pole5 can be derived from the decay widths Γ Z→ηQ+Q QX through the formulas derived in the Appendix A1 of Ref. [67], i.e., If the luminosity of a Z factory can be up to 10 35 cm −2 s −1 [43], there are about 1.4 × 10 6 η c + ccX events and 2.6 × 10 5 η b + b bX events to be produced per operation year.Moreover, the background of the Z factory is clean.Therefore, the two production processes may be studied at a high luminosity Z factory.

FIG. 4 .
FIG.4.The LO and NLO decay widths for Z → ηc +c+ c+X as functions of the renormalization scale µR.

FIG. 5 .
FIG. 5.The LO and NLO decay widths for Z → η b +b+ b+X as functions of the renormalization scale µR.

ΓFIG. 6 .
FIG.6.The LO and NLO differential decay widths dΓ/dz for Z → ηc + c + c + X.The error bars show the total uncertainties caused by the renormalization scale, the heavy quark mass and the wave function at the origin, and the total uncertainties are obtained by adding each uncertainty in quadrature.

FIG. 7 .
FIG. 7. The LO and NLO differential decay widths dΓ/dz for Z → η b + b + b + X.The error bars show the total uncertainties caused by the renormalization scale, the heavy quark mass and the wave function at the origin, and the total uncertainties are obtained by adding each uncertainty in quadrature.

FIG. 8 .FIG. 9 .
FIG.8.The LO and NLO differential decay widths dΓ/dm12 for Z → ηc + c + c + X.The error bars show the total uncertainties caused by the renormalization scale, the heavy quark mass and the wave function at the origin, and the total uncertainties are obtained by adding each uncertainty in quadrature.

FIG. 10 .FIG. 11 .
FIG.10.The LO and NLO differential decay widths dΓ/dm23 for Z → ηc + c + c + X.The error bars show the total uncertainties caused by the renormalization scale, the heavy quark mass and the wave function at the origin, and the total uncertainties are obtained by adding each uncertainty in quadrature.

TABLE I .
to the NLO level.The decay width (in unit: keV) of Z → ηc+c+c+X up to the NLO level under two different choices of µR, where the input charm quark mass is taken as the central value (i.e.mc = 1.67 GeV) and the K factor is defined as K = ΓNLO/ΓLO.

TABLE II .
The decay width (in unit: keV) of Z → η b + b + b + X up to the NLO level under two different choices of µR, where the input bottom quark mass is taken as the central value (i.e.m b = 4.78 GeV) and the K factor is defined as K = ΓNLO/ΓLO.