  We compute the ${\mathcal O}(\alpha_s^2)$ corrections to the exclusive
channels $e^+e^-\to J/\psi+\chi_{cJ}$ ($J=0,1,2$) at $\sqrt{s}=10.58$ GeV
within the nonrelativistic QCD (NRQCD) factorization framework. The validity of
NRQCD factorization at ${\cal O}(\alpha_s^2)$ has been confirmed for these
double-charmonium exclusive production processes. We analyze the impact of the
$\mathcal{O}(\alpha_s^2)$ corrections on the polarized and unpolarized cross
sections, as well as the $J/\psi$ angular distributions, which largely reduce
the renormalization scale dependence but increase the ${\mathcal O}(\alpha_s)$
NRQCD predictions to some extent for $\chi_{c0,1}$. With high numerical
accuracy, our predictions for $\sigma(J/\psi+\chi_{c1,2})$ through ${\mathcal
O}(\alpha_s^2)$ are compatible with the upper limit of the \texttt{Belle}
measurement. Although the theoretical prediction for $\sigma(J/\psi+\chi_{c0})$
is consistent with both \texttt{Belle} and \texttt{BaBar} measurements within
uncertainties, there still exists serious tension between the predicted and the
measured profiles for the $J/\psi$ angular distribution. Regarding the bright
observation prospect of the $e^+e^-\to J/\psi+\chi_{c1,2}$ channels in
\texttt{Belle 2} experiment, it is interesting to compare the future
measurements with our NRQCD predictions. The more accurate measurement of
$e^+e^-\to J/\psi+\chi_{0}$ at \texttt{Belle 2} will also help to clarify the
long-standing puzzle of $J/\psi$ angular distribution.
