Production of double heavy quarkonium at super $Z$ factory
Within the color singlet model, we calculate the exclusive production of double charmonia, double bottomonia, and double $B_c$ mesons at future super $Z$ factory. The two heavy quarkonia or $B_c$'s are either two S-wave Fock states ($^1S_0, ~^3S_1$), or one S-wave and one P-wave states ($^1P_1,~^3P_J~(J=0,1,2)$). The top three $Z^0$ propagated channels in cross sections for double charmonia are $J/ψ+h_c,η_c+χ_{c2}$, and $η_c+χ_{c0}$. For double bottomonia, they are $η_b+Υ, Υ+χ_{b2}$, and $Υ+Υ$. For double $B_c$ mesons, they are $B_c^{*+}+B_c^{*-}$,$B_c^{*+}+χ_{bc2}^-$, and $η_{bc}^++B_c^{*-}$. The cross sections of double $B_c$ 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 $σ$ as functions of CM energy $\sqrt{s}$, $σ$ as functions of the renormalization scale $μ$, the angle distributions $dσ/dcosθ$, and the $p_T$ distributions $dσ/dp_{t}$. The uncertainties come from the varying masses of constituent heavy quarks bring up to 20% corrections. We also find that the initial state radiation can bring about 30%$\sim$40% suppresions when 1\%$m_Z$ energy is losing, and cross sections can increase by about $2\sim3$ 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.