  Within the nonrelativistic QCD (NRQCD) factorization framework, we compute
the long-awaited ${\mathcal O}(\alpha_s^2)$ correction for the exclusive double
charmonium production process at $B$ factories, {\it i.e.}, $e^+e^-\to
J/\psi+\eta_c$ at $\sqrt{s}=10.58$ GeV. For the first time, we confirm that
NRQCD factorization does hold at next-to-next-to-leading-order (NNLO) for
exclusive double charmonium production. It is found that including the NNLO QCD
correction considerably reduces the renormalization scale dependence, and also
implies the reasonable perturbative convergence behavior for this process. Our
state-of-the-art prediction is consistent with the BaBar measurement within
errors.
