Production of double heavy quarkonia at super Z factory

Within the color singlet model, we calculate the exclusive production of double charmonia, double bottomonia, and double Bc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_c$$\end{document} mesons at future super Z factory. The two heavy quarkonia or Bc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_c$$\end{document}’s are either two S-wave Fock states (1S0,3S1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^1S_0, ~^3S_1$$\end{document}), or one S-wave and one P-wave states (1P1,3PJ(J=0,1,2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^1P_1,~^3P_J~(J=0,1,2)$$\end{document}). The top three Z0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z^0$$\end{document} propagated channels in cross sections for double charmonia are σ(J/ψ+hc)Z0=(454.0-59.2+39.1)×10-5fb,σ(ηc+χc2)Z0=(284.2-35.3+22.2)×10-5fb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma {(J/\psi +h_c)}_{Z^0}=(454.0^{+39.1}_{-59.2})\times 10^{-5} fb, \sigma {(\eta _c+\chi _{c2})}_{Z^0}=(284.2^{+22.2}_{-35.3})\times 10^{-5} fb$$\end{document}, and σ(ηc+χc0)Z0=(137.0-16.4+9.9)×10-5fb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma {(\eta _c+\chi _{c0})}_{Z^0}=(137.0^{+9.9}_{-16.4})\times 10^{-5} fb$$\end{document}. For double bottomonia, they are σ(ηb+Υ)Z0=(428.6-60.9+67.9)×10-4fb,σ(Υ+χb2)Z0=(230.6-22.2+23.6)×10-4fb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma {(\eta _b+\Upsilon )}_{Z^0}=(428.6^{+67.9}_{-60.9})\times 10^{-4} fb, \sigma {(\Upsilon +\chi _{b2})}_{Z^0}=(230.6^{+23.6}_{-22.2})\times 10^{-4} fb$$\end{document}, and σ(Υ+Υ)Z0=(119.9-16.6+18.4)×10-4fb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma {(\Upsilon +\Upsilon )}_{Z^0}=(119.9^{+18.4}_{-16.6})\times 10^{-4} fb$$\end{document}. For double Bc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_c$$\end{document} mesons, they are σ(Bc∗++Bc∗-)Z0=(1150-84+87)×10-3fb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma {(B_c^{*+}+B_c^{*-})}_{Z^0}=(1150^{+87}_{-84})\times 10^{-3} fb$$\end{document}, σ(Bc∗++χbc2-)Z0=(1133+22-16)×10-3fb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma {(B_c^{*+}+\chi _{bc2}^-)}_{Z^0}=(1133^{-16}_{+22})\times 10^{-3} fb$$\end{document}, and σ(ηbc++Bc∗-)Z0=(634.6-45.8+47.8)×10-3fb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma {(\eta _{bc}^++B_c^{*-})}_{Z^0}=(634.6^{+47.8}_{-45.8})\times 10^{-3} fb$$\end{document}. Here the uncertainties come from the varying masses of constituent heavy quarks, which bring up to 20% corrections. The cross sections of double Bc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_c$$\end{document} mesons are roughly one order of magnitude larger than those of the double bottomonia, and two orders of magnitude larger than those of the double charmonia. To make it helpful for experimental study, we present the total cross sections σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} as functions of CM energy s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}$$\end{document}, σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} as functions of the renormalization scale μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}, the angle distributions dσ/dcosθ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\sigma /dcos\theta $$\end{document}, and the pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_T$$\end{document} distributions dσ/dpt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\sigma /dp_{t}$$\end{document}. We also find that the initial state radiation can bring about 30–40% suppresions when 1%mZ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_Z$$\end{document} energy is losing, and cross sections can increase by about 2–3 times or decrease by an order of magnitude when adopting different potential models which becomes the major source of uncertainty. The numerical results show that it might be not optimistic for the experimental observation, but it is still far from excluded at the FCC-ee and also the CEPC running in the Z factory mode.


I. INTRODUCTION
The Large Hadron Collider (LHC) provides a direct probe of the high energy frontier of particle physics.We also have the BESIII experiment at the BEPC running at the τ -charm energy region and the Belle II experiment at the super KEKB running as a B factory.As a supplementary, we further need another electron-positron collider to make precise study on the W ± , Z 0 , and Higgs physics running at proper energy regions.The Circular Electron Positron Collider (CEPC) provides three operation modes, corresponding to the Higgs factory running at center-of-momentum (CM) energy √ s = 240 GeV, the Z factory at √ s = 91.2GeV, and the W W threshold scan around √ s ∼ 160 GeV.
For the Z-factory mode, the CEPC would have the Z 0 event yields of 7 × 10 11 [1].The first stage of Future Circular Collider (FCC) is an electron-positron collider (FCC-ee), which would operate at the Z resonance for four years.It has even larger Z 0 event yields of 7 × 10 12 , since its total integrated luminosity is about nine times as that of CEPC [2].The future electron-positron colliders would provide nice platforms for the precise study on Z 0 physics.Note that, the Electron Positron Linear Collider also has the similar GigaZ running mode [3], and a super Z factory has also been proposed by another Chinese group [4].The production of heavy quarkonium is a multiscale problem for probing the quantum chromodynamics (QCD) theory.The heavy quarkonium is a bound state consisting of a heavy quark and a heavy antiquark, i.e. the Q Q′ (Q, Q ′ = c or b quarks) bounding systems.The production of the heavy quarkonium has three different momentum scales: the heavy quark mass ∼ m Q , the momentum of the heavy quark or antiquark in the quarkonium rest frame ∼ m Q v, and the kinetic energy of the heavy quark or antiquark ∼ m Q v 2 .Here, v is the relative velocity between the constituent heavy quark and antiquark in the quarkonium rest frame.The typical values are v 2 = 0.1 for b b quarkonium, and v 2 = 0.3 for cc quarkonium.Because of m Q ≫ m Q v ≫ m Q v 2 , the creation of the constituent Q Q′ pair can be calculated perturbatively, while the hadronization of Q Q′ pair to a physical color-singlet quarkonium is nonperturbative.So the production of heavy quarkonium has been a hot topic on exploring the QCD at both the perturbative and nonperturbative energy regions, which provides rich physics on strong interation.For more information on the current status of heavy quarkonium, please refer to some reviews [5][6][7][8][9].
Considering the fact that the constituent heavy quark and antiquark inside the quarkonium have the nonrelativistic nature, we can use the nonrelativistic QCD (NRQCD) factorization framework [10,11] to factor the production of heavy quarkonium into the perturbative creation of the constituent Q Q′ pair and the nonperturbative hadronization of the heavy quark pair.The NRQCD factorization formulation consists of the short-distance coefficients and longdistance matrix elements.The short-distance coefficients stand for the hard creation of the heavy constituent Q Q pair and can be calculated perturbatively with the method of the Feynman diagrams.The long-distance matrix elements stand for the nonperturbative hadronization of the binding Q Q pair with definite J P C quantum numbers into the physical heavy quarkonium.Our calculation are in the color singlet model (CSM) [12][13][14], which require the intermediate Q Q pair and the final quarkonium have the same quantum numbers.Under this assumption, the long-distance matrix elements then can be related to the radial wave functions (or their derivatives) of the heavy quarkonium at the origin, which can be extracted from the decays of the quarkonim and can be calculated within the potential models [15][16][17][18][19][20][21][22][23][24] or potential NRQCD [25] or lattice QCD [26].Thus there are no free nonperturbative parameters in CSM.
In the calculation of the short-distance coefficients in the NRQCD framework, the squared amplitudes become complicated and lengthy especially for the production of the P -wave quarkonium.To simplify the calculation, we adopt the "improved trace technology" [27][28][29][30][31][32], which is based on the helicity amplitude method.In which, the amplitudes have been expressed directly as the linear combinations of independent Lorentz structures before the polarization sum.This method helps a lot to overcome the problem of time-consuming calculation for the production of the P -wave quarkonium.
The production of double heavy quarkonia has been studied extensively both at the LHC and B factories.The hadronic production of double J/ψ has always been a hot topic because J/ψ can be very easy to be detected by its leptonic decays and this channel can be used to explore the distribution of gluons in a proton at the LHC [33][34][35][36][37][38][39][40].In addition, the photoproduction of double J/ψ [41], the production of double heavy quarkonia through diffractive interactions [42], the production of double heavy quarkonia through photon-photon interaction [43,44], and the hadronic production of double B c mesons [45] have also been studied.The production of double heavy quarkonia through electron-positron annihilation has been widely explored at B factories.In particular, the production of J/ψ+η c at B factories once challenged the NRQCD framework.At the leading order (LO) within NRQCD formulation, its total cross section is about 2 ∼ 6 f b [46][47][48].However, the measurements at B factories by the Belle and BaBar collaborations show that the cross section is much bigger, around 20 f b [49,50].This large discrepancy between theory and experiment can be reduced by the QCD next-to-leading order (NLO) corrections [51,52] and the relativistic corrections in the NRQCD [46,53,54].Recently, this issue has been solved by the QCD next-to-next-to-leading corrections [55], which gives consistent estimate with the BaBar measurement.
The future Z factory is another platform to study the production of double heavy quarkonia through the Z 0 boson decays.Although these exclusive processes are rare, we have the advantages of the vast Z 0 events and the clear background provided by the fully reconstructed quarkonia at the electron positron collider running as a Z factory.This provide us another opportunity of studying the interplay between the perturbative and nonperturbative regimes in the hadronization of heavy quark pair into quarkonia.The production of the double charmonia at LO at such Z factory is discussed in Refs.[56,57].The QCD NLO corrections to the production of double charmonia are also studied [58,59], which shows that NLO corrections are also important near the Z 0 mass pole.The QCD NLO corrections to the paired S-wave B ( * ) c production are obtained, where the NLO corrections are small at Z 0 pole [60].The associated S-wave charmonium-bottomonium production at LO at Z 0 pole are explored in Ref. [61].In the present paper, we shall study the production of double charmonia, double bottomonia, and double B c mesons at the future Z factory, where the case of two S-wave ( 1 S 0 , 3 S 1 ) states and the case of one S-wave and one P -wave ( 1 P 1 , 3 P J (J = 0, 1, 2)) states are considered.To make our work more helpful in future experimental researches, we shall discuss the dependences of cross sections on the center-of-momentum energy and the running renormalization scale, the differential angle distribution, the transverse momentum distribution, and the uncertainties caused by the masses of constituent heavy quarks and the radial wave functions (or their derivatives) of the heavy quarkonium at the origin.
The rest of this paper is organized as follows.In Sec.II, we introduce the prescription of the production of double heavy quarkonium within the CSM framework.In Sec.III, we first present the total cross sections, then discuss their differential distributions and the uncertainties.Sec.IV is reserved for a summary.

II. FORMULATIONS
The Feynman diagrams for the exclusive production of double heavy quarkonia or B c mesons in e − (p 1 )e are displayed in Fig. 1, where H 1(2) represent the heavy quarkonia or B c 's and [n (′) ] represent the spin and orbit angular momentum states [ 1 S 0 ], [ 3 S 1 ], [ 1 P 1 ], and [ 3 P J ] (J = 0, 1, 2).Here all the γ * /Z 0 FIG.1: Feynman diagrams for electron positron annihilation into double heavy quarkonia or Bc's through the virtual photon and Z 0 boson.
intermediated Fock states are in color singlet model (CSM).The differential cross sections can be factored into the short-distance coefficients and the long-distance matrix elements [10,11], The short-distance coefficients σ describe the short-distance production of two Fock states (  ) ].
The short-distance differential cross section dσ are perturbatively calculable, where stands for the average over the spin of the initial electron and positron, and sum over the color and spin of the two Fock states.In the e − e + center-of-momentum (CM) frame, the two-body phase space can be simplified as Here, s = (p 1 +p 2 ) 2 stands for the squared CM energy.The magnitude of the 3-momentum where M 1(2) are the masses of two heavy quarkonia H 1(2) and λ[a, b, c] = (a − b − c) 2 − 4bc.The θ is the angle between the momentum p 1 of the electron and the momentum q 1 of the quarkonium The hard scattering amplitude M([n], [n ′ ]) in Eq. ( 2) can be read directly from the Feynman diagrams in Fig. 1, where j represents the number of Feynman diagrams, and l (′) are the spins of the initial leptons.The vertex V µ and the propagator D µν for the virtual photon and Z 0 propagated processes have different forms, In which, α = 1, β = 0 is for virtual photon propagated processes, and α = 0, β = 1 is for Z 0 propagated processes.e is the unit of the electric charge, e Q = 1, 2/3, −1/3 for the leptons, c quark, and b quark, respectively.g is the weak interaction coupling constant, and θ W represents the Weinberg angle.p = p 1 + p 2 is the 4-momentum of the propagator.m Z and Γ Z are the mass and the total decay width of Z 0 boson.The explicit expressions of the Dirac γ matrix chains A ν j (j = 1, 2, 3, 4) in Eq. ( 4) for the case that both two quarkonia or B c 's are the S-wave states ([ Here, S and L are the spin and the orbit angular momentum of the Fock states [n (′) ].
q are the momenta of the two constituent heavy quarks of the quarkonium H 1 [n](q 1 ) with q being the relative momentum between them and M2 q 2 − q ′ are the momenta of the two quarks of the quarkonium H 2 [n ′ ](q 2 ) with q ′ being the relative momentum between them and 2) in Eq. ( 6) have the following forms Π (S,L=0) For spin-singlet state [ 1 S 0 ], ǫ a (q k ) = 1 and γ a = γ 5 ; for spin-triplet state [ 3 S 1 ], ǫ a (q k ) = ε α (q k ) is its polarization vector and γ a = γ α is the Dirac matrix with Lorentz index α.For the S-wave states, the relative momentum q and q ′ are set to zero directly.In CSM, δ ij / √ N c is the color operator for color-singlet projector of the heavy quarkonium with N c = 3.
We shall consider the case that one of two quarkonia or ).The expressions of A ν j (j = 1, 2, 3, 4) for this case can be expressed in those for both S-wave case, In which, ε β (q 2 ) is the polarization vector of the [ 1 P 1 ] state, and ε J αβ (q 2 ) is the polarization tensor for [ 3 P J ] states with J = 0, 1, 2. The derivatives over the relative momentum q ′ β of the quarkonium or B c meson H 2 [n ′ ] will give complicated and lengthy amplitudes.The relative momentum q ′ is set to zero after taking the derivatives.
To get compact analytical expression of the complicated P -wave channels and also to improve the efficiency of numerical evaluation, we adopt the "improved trace technology" to simplify the amplitudes M([n], [n ′ ]) at the amplitude level before evaluating the polarization sum.To shorten this manuscript, we don't present the prescription.For detailed techniques and examples, one can refer to the literatures [27][28][29][30][31][32].
When manipulating the squared amplitudes 2), we need to sum over the polarization vectors of the heavy quarkonia.For the spin-triplet state [ 3 S 1 ] or the spin-singlet state [ 1 P 1 ] with 4-momentum p, the polarization sum is given by [11] Jz where J z = S z or L z for [ 3 S 1 ] and [ 1 P 1 ] states, respectively.In the case of [ 3 P J ] states, the polarization sum should be performed by the selection of appropriate total angular momentum quantum number J. The sum over polarization tensors is given by [11] The masses (units: GeV) of the constituent quarks, and the squared radial wave functions at the origin |RS(0)| 2 (units: GeV 3 ) and their first derivatives at the origin |R ′ P (0)| 2 (units: GeV 5 ) for the charmonia, bottomonia and Bc mesons within the BT-potential model [62].The uncertainties of radial wave functions (or their first derivatives) at the origin are caused by the corresponding varying quark masses.Jz

S-wave P -wave
for total angular momentum J = 0, 1, 2, respectively.At last we discuss the hadronization of the Fock states into the physical heavy quarkonia, which is described by the nonperturbative color-singlet matrix element O([n (′) ]) in Eq. ( 1).The heavy-quark spin symmetry provides relations between matrix elements for the various spin states [10], Here, v is the relative velocity between the constituent heavy quark and antiquark in the quarkonium or B c rest frame.Further, vacuum-saturation approximation together with the heavy-quark spin symmetry can be used to express all the color-singlet matrix elements in terms of the radial wave functions at the origin R S (0) for S-wave states and their first derivatives at the origin R ′ P (0) for P -wave states [10], O([ O([ Since the relativistic corrections are at order v 2 or evern higher in NRQCD, we can adopt the same values of matrix elements for both the spin-singlet and spin-triplet states.Moreover, the radial wave functions (or their derivatives) of the heavy quarkonia or B c mesons at the origin can be extracted from the decays of the quarkonia or B c 's in experiments, and be calculated within the potential models [15][16][17][18][19][20][21][22][23][24] or potential NRQCD [25] or lattice QCD [26].

III. PHENOMENOLOGY A. Input parameters
For the numerical analysis, masses of the constituent c and b quarks for heavy quarkonia and B c mesons are displayed in Table I.The mass of a heavy quarkonium is always set to be the sum of the masses of its constituent quarks , which ensures the gauge invariance of the hard scattering amplitude under the NRQCD framework.Thus we adopt different constituent quark masses for the S-and P -wave heavy quarkonia and B c mesons.In our previous work [62], the radial wave functions at the origin |R S (0)| 2 and the first derivatives of radial wave functions at the origin |R ′ P (0)| 2 for various heavy quarkonia and B c mesons have been calculated under five different potential models.In this work, we adopt the results of the Buchmüller and Tye potential model (BT-potential) [15,16,62], which are also presented in Table I.The uncertainties of radial wave functions (or their derivatives) at the origin are caused by the corresponding varying quark masses, which would be taken into consideration when we discuss the uncertainties caused by varying quark masses in Sec.III C. The LO running strong coupling constant α s is adopted, which leads to α s (2m c ) = 0.26 for charmonia and B c mesons, and α s (2m b ) = 0.18 for bottomonia.Other parameters have the following values [63]: the Fermi constant GeV, the fine structure constant α = e 2 /4π = 1/130.9,the Weinberg angle θ W = arcsin √ 0.23119, and the Z 0 boson mass m Z = 91.1876GeV and its total decay width Γ Z 0 = 2.4952 GeV.

B. Total and differential cross sections
The total cross sections for the production of double heavy quarkonia and B c mesons via e − e + → γ * /Z 0 → H 1876 GeV are listed in Table II.From this table, one could draw a conclusion that the contributions from virtual photon propagated processes are negligible in comparison with Z 0 propagated processes at future super Z 0 factory.The top three Z 0 propagated channels for the production of double charmonia, double bottomonia, and double B c mesons are In Ref. [56], Chen et.al. calculate the cross sections for the production of the double charmonia for S-and P -wave states in CSM.If the same input parameters are adopted, our estimations in Table II are consistent with theirs.
For the Z factory operation mode at CEPC, the designed integrated luminosity with two interaction point and in two years is 16 ab −1 [1].Then we can estimate the events of the production of double heavy quarkonia and B c mesons.The top three channels for double charmonia are 73 J/ψ + h c events, 45 η c + χ c2 events, and 22 η c + χ c0 events.The top three channels for double bottomonia are 686 η b + Υ events, 369 Υ + χ b2 events, and 192 Υ + Υ events.The top three channels for double B c mesons are 1.84 × 10 4 B * + c + B * − c events, 1.81 × 10 4 B * + c + χ − bc2 events, and 1.02 × 10 4 η + bc + B * − c events.We take the J/ψ + h c channel as an example to discuss its observation at CEPC.The J/ψ event can be fully reconstructed through the leptonic decay J/ψ → l + l − (l = e, µ) whose branching fraction is about 12% [63].The h c can be reconstructed through the decay chain h c → η c + γ → (hadronic decays) + γ with the (50 ± 9)% branching fraction of h c → η c + γ [63] and about total 43% branching fraction of the sixteen hadronic decays of η c1 .Thus with two years' data, we can obtain only 2 events of J/ψ + h c in the Z factory mode at the CEPC.Although the events of double bottomonia and double B c meosns are larger than the events of double charmonia by about one and two orders of magnitude respectively, the observations for the production of them are also not optimistic because of the bigger suppression factors from their decays. 2In one word, two years' operation time for Z factory mode at the CEPC is not sufficient to make the study on the production of double heavy quarkonia and B c mesons.And the observation of such decays at the CEPC with sizable sample might indicates the presence of new production mechanism, for example the contributions of color-octet Fock states might dominate.
The FCC-ee would run as a super Z factory in its first four years with the designed integrated luminosity of 150 ab −1 [2], which is about nine times of that at CEPC.Then the events of the top three channels for double charmonia, double bottomonia and double B c mesons are about 9 times of those at the CEPC.Taking the double charmonium channels as an example, we can obtain 681 J/ψ + h c events, 426 η c + χ c2 events, and 206 η c + χ c0 events.For the J/ψ + h c channel, we might reconstruct about 18 signals at the detector after considering the suppressions from the branching fractions.But if we further consider the reconstruction efficiency in experiments, the prospects are not optimistic.Anyway, the Z factory operation mode at future FCC-ee would be a better choice to study the production of double heavy quarkonia and B c mesons than at CEPC.It is worth noting that the CMS collaboration at the LHC carry out the search for Z boson decays into J/ψ and Υ pairs very recently [64].No evidence is found in the channels of Z → J/ψJ/ψ, Z → Υ(1S)Υ(1S) and Z → Υ(nS)Υ(mS), but the upper limits for the branching fractions at 95% confidence level are presented.II: J/ψ + hc, ηc + χc2, ηc + χc0, J/ψ + J/ψ, ηc + hc, and J/ψ + χc2, respectively.
In Figs.2-4, we display the total cross sections σ( √ s) as functions of the CM energy √ s, the total cross sections σ(µ) as functions of the renormalization scale µ in running strong coupling constant α s (µ), the differential angle distributions dσ/dcosθ and the p t distributions dσ/dp t at √ s = m Z the double charmonia, the double bottomonia and the double B c mesons, respectively.In those figures, we only exhibit the curves of top six Z 0 propagated processes in Table II.Note that the figure legends in those three figures are different.
In the plots for total cross sections with running CM energy √ s in Figs.2-4, one can observe the dramatic peaks of cross section caused by the Z 0 resonance clearly.In the future super Z factory, the beam energy can be controlled with a very high accuracy.The natural beam energy spread scales as σ E /E with the energy resolution σ E at the energy E. At FCC-ee, the spread is 27 MeV [65].At CEPC, the expected accuracy is the order of 100 keV for the Z factory operation [1].However, the initial state radiation (ISR) at an electron-positron collider may make the CM energy √ s deviates from the Z 0 mass pole heavily.The influence of ISR on the production rates can be identified as the deviations of √ s from m Z .In Table III, we display the top six cross sections for the double heavy quarkonia and B c mesons with the deviations of both 1%m Z and 3%m Z .It is found that the total cross sections decrease by about 30%∼40% when CM energy √ s deviates from the Z 0 pole even by only 1%m Z , which reflects the effects of ISR.The dependence of total cross sections on the renormalization scale µ in the running strong coupling α s (µ) is usually large in the leading order calculation.With its interpretation as an estimate of the contribution from missing higher order, it might be indicative.In the plots for total cross sections with running scale µ in Figs.2-4, we find that the cross sections always decrease monotonically as the scale increases with a more smooth trend at higher energy.
The differential angle distributions dσ/dcosθ for the double charmonia, double bottomonia and the double B c mesons are also displayed Figs.2-4, respectively.The distribution dσ/dcosθ can be obtained easily with the help of the differential phase space in Eq. ( 3).Here, θ is the angle between the momentum p 1 of the electron and the momentum q 1 of the heavy quarkonium.It is found that the distribution of the production of Υ + χ b0 is quite flat, while other distributions have either a slight bulge or a slight hollow.We can derive the p t distribution dσ/dp t with the help of angle distributions as follows where s is the magnitude of the 3-momentum of the heavy quarkonium.As are shown in Figs.2-4, the differential cross sections versus p t of all channels are monotonically increasing.This is understood because the p t distribution is proportional to t and the distribution dσ/dcosθ changes relatively smoothly.

C. Uncertainties
To obtain reliable estimates at leading order calculation, we discuss the main uncertainty sources of the cross sections of the double quarkonium or B c production in this subsection.For the input parameters, the Fermi constant G F , the fine structure constant α, the Weinberg angle θ W , and the mass and width of the Z 0 boson are either relatively precise or an overall factor.So we shall explore the uncertainties caused by masses of constituent heavy quarks, the radial wave functions or their first derivatives at the origin3 .We only discuss the uncertainties of the Z 0 propagated channels since the contributions of virtual photon channels are smaller by three orders of magnitude at least at super Z factory.The uncertainties of cross sections caused by both varying quark masses and the radial wave functions or their first derivatives at the origin at √ s = m Z are presented in Tables IV-VI for the double charmonia, double bottomonia, and double B c mesons, respectively.The uncertainties caused by varying masses are displayed in the second columns which are calculated in the BT-potential model [62], and the mass deviations are 0.1 GeV for m c and 0.2 GeV for m b (see explicit values in Table I).Here, the affects of the varying masses on the radial wave functions or their first derivatives at the origin are also taken into consideration, which makes them not being the overall factors any more.The values from the third to the sixth columns in Tables IV-VI are the estimates with the radial wave functions or their first derivatives at the origin calculated in other four potential models.The four models are the QCD-motivated potential with one-loop correction given by John L. Richardson (J.potential) [17], the QCD-motivated potential with two-loop correction given by K. Igi and S. Ono (I.O.potential) [18,19], the QCD-motivated potential with two-loop correction given by Yu-Qi Chen and Yu-Ping Kuang (C.K. potential) [19,20], and the QCD-motivated Coulombplus-linear potential (Cor.potential) [19,[21][22][23][24].The formula and latest values of those wave functions or their first derivatives at the origin have been discussed in our earlier work [62].It is shown that the cross sections change     In the present work, we study the production of the double heavy quarkonia or B c mesons in e − (p 1 )e + (p 2 ) → γ * /Z 0 → H 1 [n](q 1 ) + H 2 [n ′ ](q 2 ) within the CSM framework, where H 1(2) represent the heavy quarkonia or B c mesons and [n (′) ] represent the color-singlet Fock states [ 1 S 0 ], [ 3 S 1 ], [ 1 P 1 ], and [ 3 P J ] (J = 0, 1, 2).For a sound estimation at the future super Z factory, total cross sections σ as functions of CM energy √ s, σ as functions of the renormalization scale µ, the angle distributions dσ/dcosθ, and the p T distributions dσ/dp t are studied, which are shown in Figs.2-4 for the double charmonia, double bottomonia, and double B c mesons, respectively.We further discuss the uncertainties of the cross sections caused by the varying masses of constituent quarks, and the non-perturbative matrix elements under different potential models.
In Tables II, it is found that the production rates of double B c meson channels are roughly one order of magnitude greater than those of double bottomonium channels with the same Fock states, and are roughly two orders of magnitude greater than those of the double charmonium channels with the same Fock states.For CEPC running in Z factory mode, we obtain O(1) events of J/ψ + h c signal because of the big suppression factors of the fully reconstructed of both heavy charmonia.For the Z factory mode at FCC-ee, we might catch over a dozen of such rare events, but it is pessimistic for the observation in experiments after taking the reconstruction efficiency into account.However, since the color-singlet model predictions often turn out to undershoot data by one or even more orders of magnitude in many quarkonium production processes, therefore a detection of the considered processes at the FCC-ee and also the CEPC is still far from excluded.As shown in Table III, the ISR may bring about 30%∼40% suppresions when the CM energy √ s deviates from the Z 0 pole even by only 1%m Z .For the uncertainty analyses in Tables IV-VI, we find that the masses of constituent quarks can bring up to about 20% corrections to the cross sections.And the total cross sections can increase by about 2 ∼ 3 times or decrease by an order of magnitude when adopting different potential models, which would be the major source of uncertainty.

FIG. 2 :
FIG. 2: Total Cross sections versus the CM energy√ s and the renormalization scale µ, the differential angle distributions dσ/dcosθ, and the pt distributions dσ/dpt at √ s = mZ for the double charmonium production.The diamond black line, cross magenta line, dashed cyan line, solid red line, dotted blue line, and the dash-dotted green line are for the top six Z 0 propagated channels in TableII: J/ψ + hc, ηc + χc2, ηc + χc0, J/ψ + J/ψ, ηc + hc, and J/ψ + χc2, respectively.

FIG. 3 :
FIG. 3: Total cross sections versus the CM energy√ s and the renormalization scale µ, the differential angle distributions dσ/dcosθ, and the pt distributions dσ/dpt and at √ s = mZ for the double bottomonium production.The diamond black line, cross magenta line, dashed cyan line, solid red line, dotted blue line, and the dash-dotted green line are for the top six Z 0 propagated channels in TableII: η b + Υ, Υ + χ b2 , Υ + Υ, Υ + χ b1 , Υ + χ b0 , and Υ + h b , respectively.
dramatically when we choose different potential models, which bring the largest uncertainties.For the production of the double charmonia, we always obtain the maximum under the B.T. potential model and obtain the minimum under the I.O.potential model.For the production of the double bottomonia, we always obtain the maximum under the B.T. potential model and obtain the minimum under the C.K. potential model.For the production of the double B c mesons, we always obtain the maximum under the I.O.potential model and obtain the minimum under the C.K. potential model.In the Tables IV-VI, percentages in brackets are the ratios of the maximum or minimum relative to the estimates under the B.T. model.
or b quarks) with the spin and orbit angular momentum states [n] and [n ′ ], respectively.In CSM, the non-perturbative long-distance matrix elements O[n(′)] describe the hadronization of the Fock state (Q Q′ )[n(′)] into the physical heavy quarkonia H 1(2) in the same quantum states [n (′

TABLE II :
[62]s sections (units: f b) for the production of the double heavy quarkonia and Bc mesons in e + e − annihilation at √ s = mZ within the BT-potential model[62].The subscripts γ * and Z 0 are for virtual photon and Z 0 boson propagated processes, respectively.The top three channels in each column are in bold face.

TABLE III :
The deviations of cross sections (units: ×10 −3 fb) for the top six channels for the production of double charmonia, double bottomonia, and double Bc mesons when the CM energy √ s deviates from the Z 0 mass pole by 1%mZ and 3%mZ .

TABLE IV :
Uncertainties of the total cross sections (units: ×10 −5 f b) caused by constituent quark masses under B.T. model, and by four other potential models for Z 0 propagated channels of double charmonia at √ s = mZ.The percentages in brackets are the ratios of the maximum or minimum in other potential models to the estimates under the B.T. model.

TABLE V :
Uncertainties of the total cross sections (units: ×10 −4 f b) caused by constituent quark masses under B.T. model, and by four other potential models for Z 0 propagated channels of double bottomonia at √ s = mZ.The percentages in brackets are the ratios of the maximum or minimum in other potential models to the estimates under the B.T. model.

TABLE VI :
Uncertainties of the total cross sections (units: ×10 −3 f b) caused by constituent quark masses under B.T. model, and by four other potential models for Z 0 propagated channels of the double Bc mesons at √ s = mZ.The percentages in brackets are the ratios of the maximum or minimum in other potential models to the estimates under the B.T. model.