In this article, we study in detail the double-J/ψ 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 } $$
B
Z
→
J
/
ψ
+
J
/
ψ
~ 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 αs, it is observed that the prediction of $$ {\mathcal{B}}_{Z\to J/\psi +J/\psi } $$
B
Z
→
J
/
ψ
+
J
/
ψ
= ($$ {1.110}_{-0.241-0.001}^{+0.334+0.054} $$
1.110
−
0.241
−
0.001
+
0.334
+
0.054
) × 10−10, which displays a fairly steady dependence on the renormalization scale, is significantly lower than the upper limits released by CMS.