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

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/\psi$ case JXWang , there are new type Feynman diagrams, in which the $(Q\bar{Q})[^{1}S_{0}^{[1]}]$ pair is produced from two virtual gluons, need to be calculated in the $\eta_{Q}$ case. These new type Feynman diagrams do not contribute to the $J/\psi$ 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 $Z\to(Q\bar{Q})[^{1}S_{0}^{[1]}]+Q+\bar{Q}$ can be calculated through

where $M_{\rm Virtual}$ is the amplitude for the virtual correction, and $d\Phi_{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\bar{Q})[^{1}S_{0}^{[1]}]$ 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 region , 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 on-mass-shell (OS) scheme is adopted, while for the renormalization of the strong coupling constant, the modified minimal subtraction ($\overline{\rm MS}$) scheme is adopted. With this renormalization scheme, the quantities $\delta Z_{i}\equiv Z_{i}-1$ can be derived Klasen:2004tz

where $\mu_{R}$ is the renormalization scale, $\gamma_{E}$ is the Euler constant. $\beta_{0}=11C_{A}/3-4T_{F}n_{f}/3$ is the one-loop coefficient of the QCD $\beta$ function, in which $n_{f}$ is the number of active quark flavors. $\beta^{\prime}_{0}=11C_{A}/3-4T_{F}n_{lf}/3$ and $n_{lf}=3$ is the number of light-quark flavors. When $\mu_{R}\in[m_{c},m_{b})$, 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 $\mu_{R}\in[m_{b},m_{t})$, we consider the charm-quark and bottom-quark loops in the gluon self-energy but neglect the top-quark loop in the gluon self-energy, i.e. $n_{f}=n_{lf}+2=5$. For $SU(3)_{c}$ group, $C_{A}=3$, $C_{F}=4/3$ and $T_{F}=1/2$.

In the calculations, the package FeynArts feynarts is employed to generate Feynman diagrams and amplitudes, the package FeynCalc feyncalc1 ; feyncalc2 is employed to carry out the color and Dirac traces, the packages $Apart apart and FIRE fire 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 looptools numerically. The final phase-space integrations are calculated with the help of the package Vegas vegas .

The real correction to the process $Z\to(Q\bar{Q})[^{1}S_{0}^{[1]}]+Q+\bar{Q}$ comes from the process $Z(p_{0})\to(Q\bar{Q})[^{1}S_{0}^{[1]}](p_{1})+Q(p_{2})+\bar{Q}(p_{3})+g(p_{4})$. Four typical Feynman diagrams are shown in Fig.3. Compared to the $J/\psi$ case, we need to deal with new type Feynman diagrams, in which the $(Q\bar{Q})[^{1}S_{0}^{[1]}]$ pair is produced from the gluon fragmentation, such as the fourth diagram of Fig.3.

Using these Feynman diagrams, the amplitude ($M_{\rm Real}$) for the real correction can be written down directly. Then the differential decay width for the real correction can be calculated through

where $d\Phi_{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 twocutoff 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 $\delta_{s}$. Then the phase space for the real correction is divided into two regions: The soft region with $E_{4}\leq m_{{}_{Z}}\delta_{s}/2$ and the hard region with $E_{4}>m_{{}_{Z}}\delta_{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

where

and

Applying the eikonal approximation to the amplitude in the soft region twocutoff ; Bassetto:1983mvz , we obtain

Up to ${\cal O}(\delta_{s})$ corrections, the differential phase space for the soft region can be factorized as twocutoff

where $d\Phi_{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}$ Denner:1991kt ; Beenakker:2002nc , we obtain

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}}\delta_{s}/2$ for the hard region, the contribution from the hard region is finite, then $\Gamma_{\rm H}^{(Q\bar{Q})[^{1}S_{0}^{[1]}]}$ 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 $\delta_{s}$, while the sum of these two contributions should be independent to the choice of $\delta_{s}$ ($\delta_{s}$ should be small enough.). Verifying this $\delta_{s}$ independence is an important test of the correctness of the calculation. We have checked the $\delta_{s}$ independence, and have found that the results are independent of $\delta_{s}$ within the error of the numerical integration when $\delta_{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\Gamma_{Z\to\eta_{Q}+Q\bar{Q}X}$ can be obtained from $d\Gamma_{Z\to(Q\bar{Q})[^{1}S_{0}^{[1]}]+Q\bar{Q}X}$ by multiplying a factor $\langle{\cal O}^{\eta_{Q}}(^{1}S_{0}^{[1]})\rangle/\langle{\cal O}^{(Q\bar{Q})[^{1}S_{0}^{[1]}]}(^{1}S_{0}^{[1]})\rangle\approx|R_{\eta_{Q}}(0)|^{2}/(4\pi)$, where $R_{\eta_{Q}}(0)$ is $\eta_{Q}$ radial wave function at the origin, which can be calculated by using the potential model pot .

In this subsection, we give the integrated decay widths for the decay channel $Z\to\eta_{Q}+Q+\bar{Q}+X$ up to the NLO level.

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., $m_{c}=1.67\,{\rm GeV}$ and $m_{b}=4.78\,{\rm GeV}$), an analysis of the uncertainties of these decay widths will be presented later. The decay widths of $Z\to\eta_{Q}+Q+\bar{Q}+X$ up to the NLO level are given in Tables 1 and 2, where $\Gamma_{\rm NLO}$ denotes the sum of the LO contribution and the NLO corrections. The renormalization scales are set as two energy scales involved in the processes, i.e., $2m_{Q}$ and $m_{{}_{Z}}$. Tables 1 and 2 show that the NLO corrections contribute significantly to the decay widths in both cases. The NLO corrections increase the decay width of $Z\to\eta_{c}+c+\bar{c}+X$ by $\sim 34\%$ for $\mu_{R}=2m_{c}$ and $\sim 112\%$ for $\mu_{R}=m_{{}_{Z}}$; and they increase the decay width of $Z\to\eta_{b}+b+\bar{b}+X$ by $\sim 28\%$ for $\mu_{R}=2m_{b}$ and $\sim 83\%$ for $\mu_{R}=m_{{}_{Z}}$.

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\to\eta_{c}+c+\bar{c}+X(Z\to\eta_{b}+b+\bar{b}+X)$, the decay width decreases by $77\%$($57\%$) at LO, and by $63\%$($39\%$) at NLO when $\mu_{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 origin ⁴⁴4The Monte Carlo numerical integration in the calculation would lead to an error. However, the error of the numerical integration is on the order of $10^{-2}\,{\rm keV}$ for the $\eta_{c}$ case, and $10^{-3}\,{\rm keV}$ for the $\eta_{b}$ case in our calculation. Thus, the error from the numerical integration is negligible.. 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\pm 0.07\,{\rm GeV}$ and $m_{b}=4.78\pm 0.06\,{\rm GeV}$. For the radial wave functions at the origin, the authors of Ref.pot did not give an error estimate. Since the potential used in Ref.pot does not include the spin effect, the wave functions calculated in this way are accurate up to corrections of relative order $v_{Q}^{2}$. Therefore, we estimate the uncertainties by attaching an error of $30\%$ of the central value for $\eta_{c}$, and $10\%$ of the central value for $\eta_{b}$. More explicitly, we take $|R_{\eta_{c}}(0)|^{2}=0.810\pm 0.243\,{\rm GeV}^{3}$ and $|R_{\eta_{b}}(0)|^{2}=6.477\pm 0.648\,{\rm GeV}^{3}$. Then we obtain

and

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 $\eta_{c}$ and $\eta_{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.60{\rm GeV}$ ($4.72{\rm GeV}$) to $1.74{\rm GeV}$ ($4.84{\rm GeV}$), the K factor changes from $1.50$ ($1.43$) to $1.49$ ($1.44$) for the $\eta_{c}$ ($\eta_{b}$) case.

Adding the uncertainties in quadrature, we obtain

and

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 $\eta_{c}$ or $\eta_{b}$ in the processes can be defined as $z\equiv 2\,p_{0}\cdot p_{1}/m_{{}_{Z}}^{2}$. The LO and NLO differential decay widths $d\Gamma/dz$ for the processes $Z\to\eta_{c}+c+\bar{c}+X$ and $Z\to\eta_{b}+b+\bar{b}+X$ are shown in Figs. 6 and 7, respectively. Figs. 6 and 7 confirm the importance of the NLO corrections. For the $\eta_{c}$ production, the magnitude of $d\Gamma/dz$ is increased obviously at small and moderate $z$ values and decreased slightly at higher $z$ values. And for the $\eta_{b}$ production, the magnitude of $d\Gamma/dz$ is increased at all $z$ values. The uncertainties for $d\Gamma/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\Gamma/dm_{12}$ and $d\Gamma/dm_{23}$ can be measured experimentally, where $m_{12}\equiv\sqrt{(p_{1}+p_{2})^{2}}$ and $m_{23}\equiv\sqrt{(p_{2}+p_{3})^{2}}$ are invariant masses of two final-state particles. We present the LO and NLO differential decay widths $d\Gamma/dm_{12}$ and $d\Gamma/dm_{23}$ for $Z\to\eta_{c}+c+\bar{c}+X$ and $Z\to\eta_{b}+b+\bar{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\Gamma/dm_{12}$ and $d\Gamma/dm_{23}$ are changed significantly after including the NLO corrections, especially for the decay $Z\to\eta_{c}+c+\bar{c}+X$. Moreover, it is found that the NLO differential decay widths $d\Gamma/dm_{12}$ and $d\Gamma/dm_{23}$ are negative at the large $m_{12}$ or $m_{23}$ for $Z\to\eta_{c}+c+\bar{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.