  The exclusive production of double excited quarkonium is comprehensive
studied, i.e., the production of double excitedcharmonium, double excited
bottomonium, and double excited $Bc$-mesons via $e^+e^-\to \gamma^*/Z^0
\to|(Q\bar{Q'})[n]\rangle +|(Q'\bar{Q})[n']\rangle$ ($Q/Q'=c$- or $b$-quarks)
at a future $Z$ factory under the nonrelativistic quantum chromodynamics
framework, where the $[n]$~/$[n']$ represents the color-singlet heavy
quarkonium states $[^1S_0]\rangle, ~[^3S_1]\rangle, ~[^1P_1]\rangle$, and
$[^3P_J]\rangle$ ($J=0,1,2$). The "improved trace technology" is adopted for
calculating the complicated $P$-wave channels for derive the analytic
expressions at the amplitude level. According to our study, the production
rates of double heavy quarkonium are considerable at the future $Z$ factory. We
obtain the cross sections for the production of double excited charmonium for
$\sigma{(|(c\bar{c})[n]\rangle+|(c\bar{c})[n']\rangle)_{total}}=1.167^{+0.113}_{-0.164}\times
10^{-2}~fb$, the cross sections of double excited bottomonium for
$\sigma{(|(b\bar{b})[n]\rangle+|(b\bar{b})[n']\rangle)_{total}}=0.1132^{+0.0139}_{-0.0130}~fb$,
and the cross sections of double excited $Bc$-mesons for
$\sigma{(|(c\bar{b})[n]\rangle+|(b\bar{c})[n']\rangle)_{total}}=3.692^{+0.124}_{-0.097}~fb$.
The main uncertainties come from the mass of the heavy quarkonium and the
radial wave functions at the origin and their derivatives at the origin under
different potential models. The numerical results show that such the super $Z$
factory should be a good platform to study the properties of the double excited
charmonium, the double excited bottomonium , especially the double excited
$Bc$-mesons.
