Production of excited heavy quarkonia in $e^+e^- \to \gamma^*/Z^0 \to |(Q\bar{Q})[n]\rangle +\gamma$ at super $Z$ factory

Within the nonrelativistic quantum chromodynamics framework, we make a comprehensive study on the exclusive production of excited charmonium and bottomonium in $e^+e^-\to \gamma^*/Z^0 \to|(Q\bar{Q})[n]\rangle +\gamma$ ($Q=c$ or $b$ quarks) at future $Z$ factory, where the $[n]$ represents the color-singlet $n^1S_0,~n^3S_1,~n^1P_0$ and $n^3P_J$ ($n=1,2,3,4; J=0,1,2$) Fock states. The"improved trace technology"is adopted to derive the analytic expressions at the amplitude level, which is useful for calculating the complicated $nP$-wave channels. Total cross sections, differential distributions, and uncertainties are discussed in system. According to our study, production rates of heavy quarkonia of high excited Fock states are considerable at future $Z$ factory. The cross sections of charmonium for $2S$, $3S$, $4S$, $1P$, $2P$, $3P$ and $4P$-wave states are about $53.5\%$, $30.4\%$, $23.7\%$, $13.7\%$, $6.8\%$, $9.2\%$, and $9.2\%$ of that of the $1S$ state, respectively. And cross sections of bottomonium for $2S$, $3S$, $4S$, $1P$, $2P$, $3P$ and $4P$-wave states are about $39.3\%$, $12.3\%$, $14.3\%$, $7.1\%$, $3.1\%$, $2.7\%$, and $3.1\%$ of that of the $1S$ state, respectively. The main uncertainties come from the radial wave functions at the origin and their derivatives at the origin under different potential models. Then, such super $Z$ factory should be a good platform to study the properties of the high excited charmonium and bottomonium states.


I. INTRODUCTION
In comparison to the hadronic colliders like Large Hadron Collider (LHC), an electron-positron collider has some advantages, as it provides a cleaner hadronic background and the collision energy and polarization of incoming electron and positron beams can be well controlled.A super Z factory running at the energy of the Z 0 -boson mass with high luminosity L ≈ 10 34∼36 cm −2 s −1 has been proposed [1], which is similar to the GigaZ mode at an Electron-Positron Linear Collider [2] and the Circular Electron-Positron Collider (CEPC) [3].Due to the high yields of Z 0 bosons up to 7 × 10 11 at CEPC [3], it can be used for studying the production of heavy quarkonium through Z 0 decays.
The heavy quarkonium provides an ideal platform to investigate the properties of bound states, which is a multiscale problem for probing quantum chromodynamics (QCD) theory at all energy regions.Lots of data for the production of heavy quarkonium in different collisions are collected.Taking J/ψ as an example, the cross section of the inclusive production in e + e − → J/ψ + X is measured by the Bell experiment [4], the two-photon scattering in e + e − → e + e − J/ψ + X is studied by the DELPHI experiment at LEP II [5], the photoproduction in ep → J/ψ +X is explored by Zeus and H1 experiments at HERA [6,7], the hadroproduction in pp → J/ψ + X is studied by CDF experiment at Tevatron [8], and the hadroproduction in pp → J/ψ + X is widely explored by ATLAS, CMS, ALICE and LHCb experiments at LHC [9][10][11][12].Meanwhile, lots of theoretical and phenomenal efforts have been made to explain the measurements and to explore QCD.We refer the readers to some review papers to get detailed information on the status, puzzles and prospects on heavy quarkonium [13][14][15].
Considering the fact of the nonrelativistic nature of heavy quark and antiquark inside the quarkonium, the nonrelativistic QCD (NRQCD) [16,17] could be a powerful tool to study the production and decay mechanism of heavy quarkonium.In NRQCD framework, the relativistic effect with orders of v Q (v Q ≪ 1) has been separated from the nonrelativistic contributions, with v Q being the typical relative velocity between heavy quark and antiquark in the quarkonium rest frame.v 2 c ≈ 0.3 for charmonium and v 2 b ≈ 0.1 for bottomonium.Meanwhile, it divides the calculation into short-distance coefficients and the long-distance matrix elements.The short-distance coefficients describe the hard scattering of partons and can be calculated perturbatively via Feynman diagrams.The long-distance matrix elements describe the hadronization of Fock states with J P C quantum numbers into heavy quaronium and are nonperturbative parameters.
It is known that analytical expressions for the usual squared amplitudes in short-distance coefficients become complicated and lengthy for massive particles in final states especially for processes involving the P -wave Fock states.To solve the problem, the "improved trace technology" is suggested and developed [18][19][20][21], which is based on the helicity amplitudes method and deals with the trace calculation directly at the amplitude level.In this way, the amplitudes could be expressed with the linear combinations of independent Lorentz structures.In this paper, we adopt this technology to derive the analytical expression for all processes.
In previous works [22,23], the production of ground states (1S and 1P -wave) charmonium in e + e − → γ * /Z 0 → |(cc) + γ at super Z factory is studied at the leading order and next-to-leading order in strong coupling constant α s within NRQCD framework.The production of the ground states of both charmonium and bottomonium in e at Z 0 peak are explored in Ref. [24], where the contribution from initial state radiation is also considered.The production of the ground states of charmonium via virtual photon propagator in e + e − → γ * → |(cc) + γ at Bfactories are discussed in system in Refs.[25][26][27].In the present paper, we shall concentrate our attention on the production of both ground and high Fock states of both charmonium and bottomonium in e and [n 3 P J ] Fock states (n = 1, 2, 3, 4; J = 0, 1, 2).The analysis on differential distributions and the uncertainties shall be discussed.This would be a helpful support for the experimental exploration on production of those high excited heavy charmonium and bottomnium at future super Z factory or GigaZ mode at CEPC.
In the literatures [28][29][30][31], we study the production of high excited heavy quarkonium in the decay of W ± , top quark, Z 0 and Higgs boson.The numerical results show that we can obtain sizable events of heavy quarkonium of high excited [nS] and [nP ]-wave states (n ≥ 2), which implies that one can explore the special properties of those high excited states in experiments and should consider their contributions to the ground states properly.According to our study, in the processes of high excited sates could also be generated massively in comparison with the ground states.
The rest of the manuscript is organized as follows.In Section II, we introduce the calculation formalism and "new trace technology" for the processes of e + e − → γ * /Z 0 → |(Q Q)[n] + γ within the NRQCD factorization framework.In Section III, we evaluate the cross sections.The differential distributions of the cross sections and the uncertainties from various sources are studied in Sections III B and III C, respectively.The final Section IV is reserved for a summary.
The short-distance differential cross section dσ are perturbatively calculable, and the two Feynman diagrams of the processes of e − (p 1 )e + (p 2 ) → γ * /Z 0 → |(Q Q)[n] (q 1 ) + γ(q 2 ) are displayed in Fig. 1.Since the Feynman diagrams with initial state radiation can be identified in experiments, they are not considered here.The perturbative differential cross section can be expressed as where stands for the average over the spin of the initial particles and sum over the color and spin of the final particles when manipulating the squared amplitudes |M (n)| 2 .In the e − e + center-of-momentum (CM) frame, the two-body phase space can be simplified as In the second equation, we have made the integration over the δ function and the azimuth angle, and θ is the angle between the momentum p 1 of electron and the momentum q 1 of heavy quarkonium.The parameter s = (p 1 + p 2 ) 2 stands for the squared CM energy.The magnitude of the 3-dimension quarkonium momentum is , where M Q Q is the mass of heavy quarkonium.
The hard scattering amplitude M(n) in Eq. ( 2) can be read directly from the Feynman diagrams in Fig. 1.And the general form of their amplitudes can be formulated as where the index k represents the number of Feynman diagrams, and s and s ′ are the spins of the initial particles.
The vertice L µ and the propagator D µν for the virtual photon and Z 0 propagated processes have different forms, The upper and lower expressions after the big left bracket are for virtual photon and Z 0 propagated processes, respectively.In which, e is the unit of the electric charge, g is the weak interaction coupling constant, θ W represents the Weinberg angle, and m Z and Γ Z are the mass and the total decay width of Z 0 boson, respectively.
The explicit expressions of the Dirac γ matrix chains A ν k in Eq. ( 4) for the S-wave spin-singlet n 1 S 0 and spintriplet n 3 S 1 states (n = 1, 2, 3, 4) can be formulated as The first two amplitudes are for S-wave spin-singlet states and the last two are for S-wave spin-triplet states.ǫ α (q 1 ) is the polarization vector for the spin-triplet states.Π 0 q1 (q) and Π α q1 (q) are the projectors for spin-singlet states and spin-triplet states respectively, with q being the relative momentum between the two constituent quarks of heavy quarkonium.The two projectors have the following form where q 11 = q1 2 + q and q 12 = q1 2 − q are the momenta of the two constituent heavy quarks, and δ ij / √ N c is the color operator for color-singlet projector with N c = 3.For the S-wave states, the relative momentum q is set to zero directly.The vertex R ν in Eq. ( 7) is where the upper and lower expressions after the big left bracket are for the virtual photon and Z 0 propagated processes, respectively.Here e Q = 2/3 for c quark and e Q = −1/3 for b quark.We now turn to the Dirac γ matrix chains A ν k in Eq. ( 4) for the P -wave spin-singlet n 1 P 1 and spin-triplet n 3 P J states (n = 1, 2, 3, 4), which can be expressed in terms of the S-wave ones in Eq. ( 7), ; .
The first two amplitudes are for P -wave spin-singlet states and the last two are for P -wave spin-triplet states.
In which, ǫ β (q 1 ) is the polarization vector of the n 1 P 1 states and ε J αβ (q 1 ) is the polarization tensor for n 3 P J states with J = 0, 1, 2. The derivatives over the relative momentum q β in Eq. ( 10) will give complex and lengthy amplitudes.
When manipulating the squared amplitudes |M (n)| 2 , we need to sum over the polarization vectors of the heavy quarknium.For the spin-triplet n 3 S 1 states or the spinsinglet n 1 P 1 states, the polarization sum is given by [17] Jz where J z = s z or l z for n 3 S 1 and n 1 P 1 states, respectively.In the case of n 3 P J states, the polarization sum should be performed by the selection of appropriate total angular momentum quantum number.The sum over polarization tensors is given by [17] ε αβ ε (1) * for total angular momentum J = 0, 1, 2, respectively.
To get compact analytical expression of the complicated nP -wave channels and also improve the efficiency of numerical evaluation, we adopt the "improved trace technology" to simplify the amplitudes M(n) at the amplitude level before evaluating the polarization sum.To shorten this manuscript, we present its main idea below.For detailed techniques and more examples, one can refer to literatures [18][19][20][21].
Firstly, we introduce a massless spinor with negative helicity u − (k 0 ), which satisfies the following projection where k 0 is an arbitrary light-like momentum, k 2 0 = 0, and ω − = (1 − γ 5 )/2.Then we construct the massless spinor with positive helicity where k 1 is an arbitrary space-like momentum, k 2 1 = −1, and satisfies k 0 • k 1 = 0.It is easy to find that u + (k 0 ) has the projection relation where ω + = (1+γ 5 )/2.Using these two massless spinors, one can construct the massive spionrs for the fermion and antifermion, Secondly, by using the above identities, one can write down the amplitude M ±s±s ′ with four possible spin projections in the trace form directly where It is easy to check that M ±s±s ′ are orthogonal for each other.Thus, the squared amplitude can be written as where we introduce four new amplitudes Thirdly, to obtain the explicit and compact expressions as much as possible, we choose , and k µ 1 = iN 0 ε µνρσ p 1ν q 1ρ p 2σ , which leads to Then the amplitudes M i can be expressed as where e ) and The normalization factor N 0 is determined by ensuresing Thus after the three steps above, the amplitudes M i in Eq. ( 20) would be expressed by the linear combinations of some independent Lorentz structures.We finally discuss the non-perturbative matrix elements O H (n) in Eq. ( 1).They can be calculated through the lattice QCD [32], the potential NRQCD [33,34], or the potential models [29,[35][36][37][38][39][40][41].In this manuscript, we adopt the potential models to describe the non-perturbative hadronization of a (Q Q)[n] Fock state into the heavy quarkonium |(Q Q) [n] .For colorsinglet Fock states, the matrix elements are related to the Schrödinger wave function at the origin Ψ |(Q Q)[nS] (0) for the nS-wave Fock states, or the first derivative of the wave function at the origin Ψ ′ |(Q Q)[nP ] (0) for the nPwave states [16], Due to the fact that the spin-splitting effects are small, the same values of wave function for both the spin-singlet and spin-triplet Fock states are adopted in our calculation.Further, the Schrödinger wave function at the origin Ψ |Q Q)[nS] (0) and its first derivative at the origin Note that if one would take the color-octet Fock states into consideration, the color-octet NRQCD matrices are suppressed by certain orders in v Q to the corresponding color-singlet ones based on the velocity scale rules of NRQCD [13,16,42].One can also derive the values of color-octet NRQCD matrix elements by fitting the experimental measurements [43,44].

A. Input parameters
In our numerical analysis, the quark mass m Q is set to be half the mass of heavy quarkonium M Q Q/2, which ensures the gauge invariance of the hard scattering amplitude under the NRQCD framework.The masses of c and b quarks for the ground and high excited quarkonia are displayed in Table I.In our previous work [29], we calculate the radial wave functions at the origin R |(Q Q)[nS] (0) and the first derivatives of radial wave functions at the origin under five different potential models.In this work, we use the results of the Buchmüller and Tye potential model (BTpotential) [37,45], which are also presented in Table I.We will discuss the uncertainties from the radial wave functions at the origin and their derivatives at the origin under different potential models in Section III C. Note that in Table I, the uncertainties of radial wave functions at the origin and their first derivatives at the origin GeV 3 ) and their first derivatives at the origin ) within the BT-potential model [29].Uncertainties of radial wave functions at the origin and their first derivatives at the origin are caused by the corresponding varying quark masses.
GeV under the BT-potential model.Percentages in brackets are ratios relative to the ground state.are caused by the corresponding varying quark masses.It tells us that the evaluation of cross sections of high excited Fock states (n = 2, 3, 4) are more than simply the non-perturbative matrix elements in the calculation for the ground state (n = 1).The nonperturbative matrix elements depend on the heavy quark masses.Other parameters have the following values [46]: the mass of Z 0 boson m Z = 91.1876GeV and its total decay width Γ Z 0 = 2.4952 GeV, the Fermi constant GeV, the Weinberg angle θ W = arcsin √ 0.23119, and the fine structure constant α = e 2 /4π = 1/130.9.
GeV under the BTpotential model.Percentages in brackets are ratios relative to the ground state.
The total cross sections for the production of heavy quarkonia via e 1876 GeV are listed in Tables II and III for virtual photon γ * and Z 0 propagated processes, respectively.The percentages in brackets are ratios of high excited states (n = 2, 3, 4) relative to the ground state (n = 1).Here we adopt the BT-potential model to evaluate the non-perturbative hadronic matrix elements [29].It is worth noting that, there are no estimations on the σ ) via the virtual photon propagated processes in Table II because they break up the conservation of C parity.In Refs.[22,23], Chen et.al. calculate the cross sections for 1S and 1P -wave charmonium in e − e + → γ * /Z 0 → |(cc)[n] + γ at leading and next-to-leading order accuracy in strong coupling constant α s .If the same input parameters are adopted, our estimations are consistent with theirs at leading order.
Since the units in Table III are two orders larger than units in Table II, the contributions from the virtual photon processes are negligible at future super Z factory.In Table III for Z 0 propagated processes, it is found that where Q = c or b quarks.For bottomonium |b b[n] , the cross sections of n 1 P 1 Fock state for all n = 1, 2, 3, 4 are quite close to those of the n 3 P 1 Fock state at the same nth level.It is worth noting that in Ref. [24], they considered the contribution from initial state radiation and found that σ(|(b b[ Let's take a closer look at the cross sections of the high excited states in Table III.When using [nS] to represent the sum of cross sections of n 1 S 0 and n 3 S 1 , and [nP ] to represent the sum of cross sections of n 1 P 1 and n 3 P J (J = 0, 1, 2) at the same nth level, we have Then at the future Z factory or CEPC in GigaZ mode running at CM energy √ s = m Z with high luminosity, we can obtain sizable events to study both ground and high excited heavy quarkonia.We can obtain the events in one operation year simply by multiplying the cross sections in Tables II and III by the luminosity L ≈ 10 36 cm −2 s −1 ≈ 10 4 f b −1 year −1 .In Figs. 2 and 3, we display the total cross sections versus the CM energy √ s for ground states |(cc) [1] and |(b b) [1] respectively, where [1] stands for 1 1 S 0 , 1 3 S 1 , 1 1 P 1 , and 1 3 P J -wave states (J = 0, 1, 2).They show explicitly the contributions of γ * and Z 0 propagated processes from √ s = 10 GeV to 140 GeV.Around the Z 0 peak, the Z 0 propagated processes dominate without any doubts.The curves of total cross sections versus In Fig. 4, differential distributions dσ/dcosθ for ground states |(cc) [1] and |(b b) [1] are displayed, where [1] stands for 1 1 S 0 , 1 3 S 1 , 1 1 P 1 , and 1 3 P J -wave states (J = 0, 1, 2).Here, θ is the angle between the momentum p 1 of electron and the momentum q 1 of the heavy quarkonium.It is shown that the Z 0 propagated processes and the corresponding virtual photon propagated ones have similar line shapes.We also find that dσ/dcosθ approaches its maximum when the heavy quarkonium and the electron running in the same direction or backto-back for both S-wave and P -wave states.The curves of differential cross sections dσ/dcosθ for high excited states |(cc) [ (left), the Z 0 boson (middle) and the sum of previous two (right).The diamond black line, cross magenta line, dashed cyan line, solid red line, dotted blue line, and the dash-dotted green line are for distribution dσ/dcosθ is set to be which can be easily obtained with the differential phase space of Eq. ( 3), then the distribution dσ/dp t can be obtained by where is the magnitude of the momentum of the heavy quarkonium.We present the transverse momentum p t distributions for the cross sections in Fig. 5 for ground states |(cc) [1] and t and values of the function f (cosθ) changes smoothly, dσ/dp t shall increase with the increment of transverse momentum p t .The curves of differential cross sections dσ/dp t for high excited states TABLE IV: Uncertainties of total cross sections (units: ×10 −4 fb) caused by varying the masses as shown in Table I for γ * propagated processes.Note, effects of uncertainties of radial wave functions at the origin and their first derivatives at the origin caused by varying masses are also considered.

C. Uncertainty analysis
For the leading-order calculation, the main uncertainty sources of cross sections include the Fermi constant G F , the Weinberg angle θ W , the fine-structure constant α, the mass and width of the Z 0 boson, the masses of constituent quarks, and the non-perturbative matrix elements.Since parameters G F , θ W , α and the mass and width of the Z 0 boson are either an overall factor or an relatively precise value, we will not discuss uncertainties caused by them.In this subsection, we will explore uncertainties caused by masses of constituent quarks, the non-perturbative matrix elements, and deviation of CM energy √ s away from m Z .
The uncertainties of cross sections caused by varying the masses of constituent quarks by 0.1 GeV for m c and 0.2 GeV for m b (as shown in Table I ) at the CM energy √ s = 91.1876GeV are presented in Tables IV and V for virtual photon γ * and Z 0 propagated processes, respectively.It worth noting that the effects of uncertainties of radial wave functions at the origin and their first derivatives at the origin caused by varying masses are also taken into consideration.It is found that the wave functions at the origin and their derivatives at the origin increase as quark masses increase.But, we find that the short-distance coefficients decrease along with the increasement of quark masses.The overall effect is that the cross sections decrease with the increment of the quark masses.
We adopt four other potential models to estimate the uncertainties caused by the wave functions at the origin and their first derivatives at the origin in Tables VI and VII for Z 0 propagated processes for charmonium and bottomonium, respectively.The four models are QCDmotivated potential with one-loop correction given by John L. Richardson (J.potential) [47], QCD-motivated potential with two-loop correction given by K. Igi and S. Ono (I.O.potential) [48,49], QCD-motivated potential with two-loop correction given by Yu-Qi Chen and Yu-Ping Kuang (C.K. potential) [40,49], and the QCDmotivated Coulomb-plus-linear potential (Cor.potential) [35,36,[49][50][51].The formula and latest values of those wave functions at the origin and their first derivatives at the origin can be found in our earlier work [29].In Tables VI and VII, the contributions from four P-wave states (n 1 P 1 , n 3 P J with J = 0, 1, 2) are summed up.It is shown that the cross sections change dramatically when we choose different potential models.I for Z 0 propagated processes.Note, effects of uncertainties of radial wave functions at the origin and their first derivatives at the origin caused by varying masses are also considered.
models.In Tables VI and VII, percentages in brackets are the ratios of the minimum or maximum relative to the estimates under the B.T. model.For the uncertainties of total cross sections caused by the deviation of CM energy √ s away from m Z , one can have a visual impression in Figs. 2 and 3.It is shown that the cross sections decreases dramatically with the deviation of CM energy √ s away from m Z .To obtain a quantitative impression, we display the uncertainties caused by the deviation of CM energy √ s away from m Z by 1% and 3% for the Z 0 propagated process with n = 1 in Table VIII.
and [n 3 P J ] Fock states (n = 1, 2, 3, 4; J = 0, 1, 2).The"improved trace technology", which disposes the Dirac matrices at the amplitude level, is helpful for deriving compact analytical results especially for the complicated P -wave processes with massive spinors.The total cross sections σ( √ s) and differential distributions dσ/dcosθ and dσ/dp t for all n = 1 Fock states are studied in detail.For a sound estimation, we further study the uncertainties of the cross sections caused by the varying mass of c and b quarks, the non-perturbative matrix elements under five potential models, and deviation of CM energy √ s away from m Z .In addition to the ground states, it is found that the production rates of high excited Fock states of charmonium and bottomonium are considerable in the processes of e + e − → γ * /Z 0 → |(Q Q)[n] + γ at super Z factory with high luminosity L ≈ 10 36 cm −2 s −1 .The cross sections of charmonium for 2S, 3S, 4S, 1P , 2P , 3P and 4P -wave states are about 53.5%, 30.4%, 23.7%, 13.7%, 6.8%, 9.2%, and 9.2% of that of the 1S state, respectively.And cross sections of bottomonium for 2S, 3S, 4S, 1P , 2P , 3P and 4P -wave states are about 39.3%, 12.3%, 14.3%, 7.1%, 3.1%, 2.7%, and 3.1% of that of the 1S state, respectively.Then, such a super Z factory could provide a useful platform to study the high excited charmonium and bottomonium.In addition, we find that cross sections change dramatically when adopting different potential models, which would be the major source of uncertainty.And the deviation of CM energy √ s away from Z 0 pole at future super Z factory will also have great influence on the production rates.
For the production of |(cc)[n] in Table VI, we always obtain the minimum under the I.O.potential model, and obtain the maximum under the B.T. potential or J. potential models.While for the production of |(b b)[n] in IV. CONCLUSIONSIn the present work, we make a comprehensive study on the high excited states of the |(cc)[n] and |(b b)[n] quarkonium production in e + e − → γ * /Z 0 → |(Q Q)[n] + γ within the NRQCD factorization framework at future Z factory, where [n] stands for [n 1

TABLE I :
Masses (units: GeV) of the constituent quark and radial wave functions at the origin |R 2 , 16.12 +1.28

TABLE II :
Cross sections (units:
Table VII, we obtain the minimum under C.K. or I.O.potential models, and obtain the maximum under the B.T., or Cor.potential

TABLE V :
Uncertainties of total cross sections (units: ×10 −2 fb) caused by varying the masses as shown in Table

TABLE VII :
Uncertainties of total cross sections (units: ×10 −2 fb) caused by five different potential models for |(b b)[n] quarkonium in e + e − → Z 0 → |(b b)[n] + γ.Percentages in brackets are the ratios of the minimum or maximum relative to the estimates under B.T. model.