Comprehensive studies of $\Upsilon$ inclusive production in $Z$ boson decay

In this paper, we present a comprehensive study of $\Upsilon$ inclusive production in $Z$ boson decay, including the first complete next-to-leading-order calculations of the color-octet (CO) contributions. With the inclusion of the newly-calculated remarkable QCD corrections, the CO processes exhibit crucially phenomenological influence on the existing predictions built on the color-singlet mechanism. We also include the exhaustive evaluations of the feed-down contributions, which remained ignored in the literature, and find them to be considerable. Summing up all the contributions, the $\mathcal{B}_{Z \to \Upsilon(nS)+X}$ still notably undershoot the data released by the L3 Collaboration.


Introduction
Heavy-quarkonium production in Z boson decay, which has attracted much attention on both theoretical and experimental sides in the past decades , can offer independent test of the quarkonium-production mechanism and provide references to distinguish different models. A well-known example is that the measurements of J/ψ inclusive production in Z decay apparently exceed the predictions given by the color-singlet (CS) mechanism but coincide with that built on the non-relativistic QCD (NRQCD) framework [8][9][10]34], which is regarded as a solid evidence to favor the NRQCD factorization.
In addition to inclusive J/ψ production, L3 Collaboration also measured the total decay width of Z → Υ(nS) + X [16], (1.1) The leading-order (LO) calculations of the CS process Z → bb[ 3 S [1] 1 ]+b+b give a prediction B Z→Υ(1S+2S+3S)+X ∼ 10 −6 [3,4]. Subsequently, Li et al. [19] evaluated the next-toleading-order (NLO) QCD corrections to this process, and pointed out that the higherorder terms in α s could give rise to a 24 − 100% enhancement, depending on the choice of the renormalization scale (µ r ). Recently, we accomplished the NLO calculations of Z → bb[ 3 S [1] 1 ] + g + g, which is of the same order in α s as Z → bb[ 3 S [1] 1 ] + b +b, discovering that this process is indispensable [33]. Summing over all the above contributions, the CS prediction of B Z→Υ(1S+2S+3S)+X is just about 10 −5 and therefore incompatible with the L3 measurements.
In order to fill the huge gap between theory and experiment, it is urgent to evaluate the CO contributions, among which only the result for the Υ( 3 S [8] 1 ) production via gluon fragmentation 1 , Z → q +q + g * with g * → bb[ 3 S [8] 1 ], has been given [9]. The author of Reference [9] claimed that the decay width of Z → Υ + X can be enhanced by several times with the inclusion of the above gluon-fragmentation contribution alone, which, however, can be attributed to the employment of a very large value of O Υ ( 3 S [8] 1 ) . In some recent references [35][36][37], this value has been renewed, and found to be more than one order of magnitude smaller than that used in Reference [9]. Further, we will find in this paper that the gluon-fragmentation part actually does not dominate the whole CO contributions, as a result, a complete calculation of the CO processes is necessary. To make a comprehensive study, the feed-down contributions from higher excited states, namely Υ(2S, 3S) and χ b (1P, 2P, 3P ), need also to be counted. As a matter of fact, this part is nonnegligible. Taken together, in this paper we will revisit the Υ inclusive production in Z decay by including the first complete NLO studies of the CO processes and the exhaustive evaluations of the feed-down contributions.
Note that, in comparison with J/ψ, Υ may be even more suitable for the investigations of heavy-quarkonium inclusive production in Z decay. The large b-quark mass will make the Υ decay products more energetic and thereby more easily detectable. What is more important is, in Z → J/ψ + X, the significant (even dominant) b-hadron feed-down contributions will impose great obstacles to the extraction of prompt J/ψ production rates; that is to say, the value of B b→J/ψ+X severely affects the precision of fitting B Z→J/ψprompt+X [16]. However, there is no such b-hadron feed down contributing in Z → Υ + X, which would be rather beneficial to achieve a precise measurement of B Z→Υprompt+X . On the theoretical side, the large mass of b quark generally results in a better convergent perturbative series over the expansion of α s and v 2 . Considering the large uncertainty of B Z→Υ(1S+2S+3S)+X and the lack of definitely measured B Z→Υ(nS)+X as listed in equation (1.1), which are primarily attributed to the low Υ production rates at LEP, it is crucial to reperform the measurements with better precision at colliders equipped with much higher luminosities, such as LHC or some planned Z factories. Our state-of-the-art predictions would pave the way for comparisons with the future measurements.
The rest of the paper is organized as follows: Section 2 is an outline of the calculation formalism. Then, the phenomenological results and discussions are presented in Section 3. Section 4 is a concluding remark. 1 Actually, the process Z → q+q+g * ; g * → bb[ 3 S [8] 1 ] is only one part of the gluon-fragmentation processes; the other one of the same order in αs is the loop-induced process Z → g + g * ; g * → bb[ 3 S [8] 1 ], which was uncalculated in the literature and which will be computed in our present NLO calculations.

Calculation formalism
In the context of NRQCD factorization [34,38], the decay width of Z → Υ(χ b ) + X 2 can be factorized as J ; for the χ b case, up to O(v 2 ), n can take 3 P

CS framework
(2.2) NLO QCD corrections to the first one have been carried out by Li et al. [19], and thus we straightforwardly employ their results. Regarding the gluon-radiation process in equation 1 ] + g + g, we reperform the NLO calculations by taking the presently used parameters, based on the formalism described in our recent paper [33].
For n = 3 P

[1]
J , the LO processes follow as (2.3) From Reference [31] we learn that the χ b feed-down contributions to Υ production through the two CS processes in equation ( 1 ) 3 , which are beyond the scope of this article. Combining the two points, we compute the two 3 P [1] J processes in equation (2.3) only at the tree-level accuracy, which would not contaminate the precision of our NRQCD predictions of Γ Z→Υ+X .

CO framework
In the case of n = 1 S [8] 0 , 3 S 2 Z → χ b + X should be also calculated so as to obtain the χ b feed-down contributions. 3 The tree-level process Z → bb[ 3 P [1] J ] + g + g by itself has contained soft singularities, which arise from the attachment of a soft gluon to the P-wave quarkonium and which can only be canceled by introducing Γ Born refers to the LO processes;Γ Virtual andΓ Real represent the virtual and real corrections, respectively.Γ NLO * denotes the heavy quark-antiquark pair associated processes, i.e., Z → bb[ 1 S , which are free of divergences and thereby can be directly computed using the standard Monte-Carlo integration techniques. In the following, we just briefly describe the formalism for calculating the virtual and real corrections.
1 ) should be additionally added to eliminate the residually soft singularities.

Virtual corrections
The virtual corrections consist of the one-loop contributions (Γ Loop ) and the counterterm contributions (Γ CT ), whose typical diagrams are shown in figure 1;Γ Virtual can then be expressed asΓ We utilize the dimensional regularization with D = 4 − 2 to isolate the ultraviolet (UV) and infrared (IR) divergences. The on-mass-shell (OS) scheme is employed to set the renormalization constants for the heavy quark mass (Z m ), heavy quark filed (Z 2 ), and gluon filed (Z 3 ). The modified minimal-subtraction (M S) scheme is used for the QCD gauge coupling (Z g ). The renormalization constants are taken as where γ E is the Euler's constant, is an overall factor in our calculation, β 0 = 11 3 C A − 4 3 T F n f is the one-loop coefficient of the β function, and β 0 = 11 3 C A − 4 3 T F n f l . n f (= 5) and n f l (= n f − 2) are the numbers of active-quark flavors (including c and b quarks as in Reference [39]) and light-quark flavors, respectively. In SU(3), the color factors are given by T F = 1 2 , C F = 4 3 , and C A = 3. To calculate the D-dimension trace of fermion loops involving γ 5 , under the scheme described in [19,22,40,41], we choose the same starting point (Z-vertex) to write down all the amplitudes without implementation of cyclicity.

Real corrections
The real corrections to Z → bb[n]+g (n = 1 S as illustrated in figure 2. "q" and "u g " represent the light quarks (u, d, s) and ghost particles, respectively. The phase-space integration of the three processes in equation (2.8) will generate IR singularities, which can be isolated by slicing the phase space into different regions, namely the two-cutoff slicing strategy [42]. By introducing two small cut-off parameters (δ s and δ c ) to decompose the phase space into three parts,Γ real can then be written asΓ Γ S is the soft term arising only from the process of Z → bb[n] + g + g;Γ HC denotes the hardcollinear term, which originates from all the three processes in equation (2.8 (2.10) Then one can obtain u is finite, and its analytical expression can be found in Reference [43].
0 ) are adopted. 5 Actually, in our NLO calculations, only the axial-vector part of Z → bb[ 3 P [8] J ] + g + g comprises the so-called residually soft singularities.  Figure 3: Verifications of the independence on the soft (δ s ) and collinear (δ c ) cut-off parameters.
In our calculations, we use our M athematica-F ortran package with the implementation of FeynArts [44], FeynCalc [45], FIRE [46], and Apart [47]. This package has been employed to preform QCD corrections to several heavy-quarkonium related processes [33,48,49]. Simultaneously, we apply another independent package, Feynman Diagram Calculation (F DC) [50], to compute all involved processes, and acquire the same numerical results. Independence on the cut-off parameters has been scrutinized, as is shown in figure 3 6 .
As a crosscheck, by taking the same input parameters, we have reproduced the NLO results of σ(e + e − → cc[ 1 S [8] 0 , 3 P

Input parameters
The input parameters entering our calculations are taken as As to the NRQCD LDMEs, we utilize three typical sets of these parameters [35][36][37] 6 The process Z → q +q + g * ; g * → bb[ 3 S [8] 1 ], which is divergence free and therefore does not depend on δs,c, contributes significantly toΓ HC ( 3 S With a view to applications of M 0 and M 1 , we set O Υ(nS) ( 1 S [8] 0 ) to ζM Υ(nS) 0,r 0 , and correspondingly 0,r 0 , and vary ζ between 0 and 1 around the default value 1/2.
J ) , are related to the radial wave functions and their first derives at the origin (n, m = 1, 2, 3) by the following formula:

SDCs
Before going further, we first take a look at the SDCs. To begin with, we summarize the CS SDCs in figure 4. Inspecting the above two figures, the QCD corrections to the two 3 S ; however, further multiplying by 3 4π |R χ bJ (mP ) (0)| 2 and B χ bJ →Υ+γ would eventually result in an insignificant χ b feed-down contribution to Υ productions. That is to say, the 3 S The CO SDCs with respect to µ r are drawn in figure 5. One can observe that the NLO QCD corrections to Z → bb[ 1 S [8] 0 , 3 P [8] J ] + X could enhance their LO results by about a steady factor of 2. The NLO * contributions given by Z → bb[ 1 S [8] 0 , 3 P  Figure 4: SDCs of the CS processes. The subscript "gg" stands for Z → bb[ 3 S [1] 1 ] + g + g [33], and "bb" for Z → bb[ 3 S [1] 1 ] + b +b [19], withΓ gg /Γ bb elucidating the ratio of these two channels. In derivingΓ 3 P 1 state, the high-order terms in α s can magnify its LO results to a extremely large extent, by about 25-30 times, which can be mainly ascribed to the gluon-fragmentation diagrams that first appear at the NLO level. By the same token, the 3 S 1 NLO * processes further enlarge its NLO results by about 50%. The rightmost figure in the second column suggests that the SDC of 3 P [8] J is comparable to that of 3 S can reach up to about 30%, increasing towards higher µ r . AlthoughΓ 1 S [8] 0 is only about 10% in magnitude ofΓ 3 S [8] 1 , the huge value of O Υ ( 1 S [8] 0 ) [35][36][37] would compensate for its smallness and provide a even larger contribution than that of 3 S [8] 1 . Therefore, the existing theoretical studies, which just concentrate on the process of Z → q +q + g * ; g * → bb[ 3 S [8] 1 ], are indeed far insufficient to provide us with a thorough CO result. Table 1: Decay widths of Z → Υ(3S) + X (in units of KeV). The subscripts "Dr" and "Fd" refer to the direct-production processes and feed-down effects, respectively. m b = 4.75 GeV and µ r = m Z /2.  Table 2: Decay widths of Z → Υ(2S) + X (in units of KeV). The subscripts "Dr" and "Fd" refer to the direct-production processes and feed-down effects, respectively. m b = 4.75 GeV and µ r = m Z /2.  Table 3: Decay widths of Z → Υ(1S) + X (in units of KeV). The subscripts "Dr" and "Fd" refer to the direct-production processes and feed-down effects, respectively. m b = 4.75 GeV and µ r = m Z /2. The NRQCD predictions of Γ Z→Υ(nS)+X are presented in tables. 1, 2, and 3, and in figure 6. We find 1) The CO contributions exhibit crucial influence on the theoretical predictions, the amount of which relying on the LDMEs choices. Taking the LDMEs of Refs. [36] and   Figure 6: Decay widths of Z → Υ(nS) + X as a function of µ r , obtained by the LDMEs of Refs. [35][36][37]. The subscript "NR" means the sum of the CS and CO contributions. "Dr" and "Fd" refer to the direct-production processes and feed-down effects, respectively. m b = 4.75 GeV and µ r = m Z /2.

Decay widths
[37], inclusion of the CO processes can increase the CS predictions (Γ Dr CS ) by about 20 − 30%. When the LDMEs of Reference [35] is adopted, Γ Dr CO can magnify Γ 3S Dr CS by about 40%; however, for Υ(2S) and Υ(1S) production, this kind of enhancements are mild.
Taken together, our newly-calculated CO and feed-down contributions would introduce a 30 − 45% enhancement to the existing theoretical results, manifesting the necessities of our calculations. By comparing with the L3 measurements, one can observe that the newest NRQCD predictions appear to be still at variance with the L3 data. To be specific, the predicted B Z→Υ(1S)+X is about 5 times smaller in magnitude than the experimental upper limit, and this value can rise up to about 30 for Υ(2S, 3S); in the case of B Z→Υ(1S+2S+3S)+X , the central value of the L3 measurement is about 7 times bigger in amount than our predictions. It is worth noting that the L3 Collaboration used only 6 reconstructed events to fit B Z→Υ(1S+2S+3S)+X , which is responsible, in part, for its large measuring uncertainties; moreover, the definite values of B Z→Υ(1,2,3S)+X , rather than the upper limits, have not yet been measured. In this sense, perhaps it is still premature to draw a decisive conclusion concerning the consistency or inconsistency of the NRQCD predictions with the L3 data, requiring future measurements at a collider with much higher luminosity, such as LHC or some planned Z factories. At last, we analyze the uncertainties of the predictions due to the choices of the renormalization scale µ r , the CO LDMEs, and the b-quark mass m b .

Summary
In order to tackle the great discrepancies between the theoretical results and the L3 data of Γ Z→Υ+X , we in this manuscript revisit the inclusive Υ production in Z decay by including the complete evaluations of the CO contributions at the QCD NLO accuracy for the first time. Our results show that the newly-calculated QCD corrections to the CO LO processes are significant, subsequently resulting in a large enhancement to the existing predictions given by the CS model. Besides, the feed-down contributions are considered exhaustively, and found to considerably enhance the decay width. Taking into account all theses contributions, the discrepancies between theory and L3 data remain conspicuous.