  In this article, we study in detail the double-$J/\psi$ yield through $Z$
decay at the next-to-leading-order (NLO) QCD accuracy within the
nonrelativistic QCD factorization. At the tree level, the pure QCD diagrams
predict a branching ratio of $\mathcal{B}_{Z \to J/\psi+J/\psi} \sim 10^{-12}$;
however, the inclusion of the QED diagrams would augment this prediction by
approximately 2-3 orders of magnitude. After incorporating the QCD corrections,
the QCD results exhibit a considerable increase, whereas the QED results
undergo a substantial reduction. Combining the QCD and QED contributions at NLO
in $\alpha_s$, it is observed that the prediction of $\mathcal{B}_{Z \to
J/\psi+J/\psi}=(1.110^{+0.334+0.054}_{-0.241-0.001})\times 10^{-10}$, which
displays a fairly steady dependence on the renormalization scale, is
significantly lower than the upper limits released by CMS.
