${\mathcal O}(\alpha_s^2)$ corrections to $J/\psi+\chi_{c0,1,2}$ production at $B$ factories

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.