Within the nonrelativistic QCD framework, we calculate the double heavy quarkonium (double charmonium and double bottomonium) production at a future super $Z$ factory and at the Circular Electron Positron Collider (CEPC) and Electron-Positron Future Circular Collider (FCC-ee). The complete color-octet (CO) channels in the ${\gamma }^{*}/Z^{*}$-propagated process are considered along with the color-singlet (CS) channels at leading order of ${\alpha }_s$. We found that the contributions of CO states to the total cross section are significant or dominant for many processes within the energy region at a $Z$ factory and at the CEPC and FCC-ee. The experimental measurements will help us to verify the CO mechanism. Among these CO channels, the gluon fragmentation process (into ${{}^{3}S}_1^8$ states at leading order) is most important. The comparison between the theoretical results and future data may give constraints to the CO matrix elements. Additionally, we consider the relativistic corrections to both the CS and CO channels which decrease the cross sections significantly. Specially, the $K$ factors are about 0.5 for most charmonium channels. We get estimates of the events for double heavy quarkonium production. The final events of $J/\psi +{\eta }_c$, $J/\psi +J/\psi$, and ${\eta }_c+{\eta }_c$ production, with $J/\psi$ reconstructed by lepton pair and ${\eta }_c$ reconstructed by hadronic decays, would be (1, 8, 4) and (11, 77, 37) at the CEPC (2-yr) and at the FCC-ee (4-yr) for the $Z$ factory mode, respectively.