  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/\psi+h_c,\eta_c+\chi_{c2}$, and
$\eta_c+\chi_{c0}$. For double bottomonia, they are $\eta_b+\Upsilon,
\Upsilon+\chi_{b2}$, and $\Upsilon+\Upsilon$. For double $B_c$ mesons, they are
$B_c^{*+}+B_c^{*-}$,$B_c^{*+}+\chi_{bc2}^-$, and $\eta_{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 $\sigma$ as functions of CM energy $\sqrt{s}$,
$\sigma$ as functions of the renormalization scale $\mu$, the angle
distributions $d\sigma/dcos\theta$, and the $p_T$ distributions
$d\sigma/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.
