  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.
