  The failure of observing the $e^+ e^- \to J/\psi+J/\psi$ events at $B$
factories to date is often attributed to the significant negative
order-$\alpha_s$ correction. In this work we compute the ${\cal O}(\alpha^2_s)$
correction to this process for the first time. The magnitude of the
next-to-next-to-leading order (NNLO) perturbative correction is substantially
negative so that the standard NRQCD prediction would suffer from an unphysical,
negative cross section. This dilemma may be traced in the fact that the bulk
contribution of the fixed-order radiative corrections stems from the
perturbative corrections to the $J/\psi$ decay constant. We thus implement an
improved NRQCD factorization framework, by decomposing the amplitude into the
photon-fragmentation piece and the non-fragmentation piece. With the measured
$J/\psi$ decay constant as input, which amounts to resumming a specific class
of radiative and relativistic corrections to all orders, the
fragmentation-induced production rate can be predicted accurately and serves a
benchmark prediction. The non-fragmentation type of the amplitude is then
computed through NNLO in $\alpha_s$ and at lowest order in velocity. Both the
${\cal O}(\alpha_s)$ and ${\cal O}(\alpha^2_s)$ corrections in the interference
term become positive and exhibit a decent convergence behavior. Our finest
prediction is $\sigma(e^+ e^- \to J/\psi+J/\psi)= 2.13^{+0.30}_{-0.06}$ fb at
$\sqrt{s}=10.58$ GeV. With the projected integrated luminosity of 50 ${\rm
ab}^{-1}$, the prospect to observe this exclusive process at \texttt{Belle} 2
experiment appears to be bright.
