  In this paper, we employ the nonrelativistic QCD factorization to conduct a
comprehensive examination of the $Z$ boson decay into a pair of $\Upsilon$
mesons, achieving accuracy at the next-to-leading-order (NLO) in $\alpha_s$.
Our calculations demonstrate that the QED diagrams are indispensable in
comparison to the pure QCD diagrams, and the implementation of QCD corrections
markedly enhance the QCD results, whereas it substantially diminish the QED
results. To ensure consistency with the experimental methodology, we have taken
into account the feed-down transitions originating from higher excited states,
which exhibit significant relevance. Combining all the contributions, we arrive
at the NLO prediction of $\mathcal{B}_{Z \to \Upsilon(nS)+\Upsilon(mS)} \sim
10^{-11}$, which is notably lower than the upper limits set by CMS.
