  Within the nonrelativistic QCD (NRQCD) factorization framework, we have
computed the ${\mathcal O}(\alpha_s^2)$ corrections to the exclusive production
of $P$-wave spin-triplet charmonia $\chi_{cJ}(J=0,1,2)$ accompanied by a hard
photon at $B$ factories. For the first time, we have explicitly verified the
validity of NRQCD factorization for exclusive $P$-wave quarkonium production at
two-loop order. Unlike the $\chi_{cJ}$ electromagnetic decay processes, the
$\mathcal{O}(\alpha_s^2)$ corrections are found to be smaller than the
$\mathcal{O}(\alpha_s)$ corrections in all three channels $e^+e^-\to
\chi_{c0,1,2}+\gamma$. In particular, the ${\mathcal O}(\alpha_s^2)$
corrections appear moderate for $\chi_{c1}$ case, and insignificant for
$\chi_{c0}$. In addition, the next-to-next-to-leading order (NNLO) predictions
for the production rates of $\chi_{c0,1}+\gamma$ are insensitive to the
renormalization and factorization scales. All of these features may indicate
that perturbative expansion in these two channels exhibits a decent convergence
behavior. By contrast, both the ${\mathcal O}(\alpha_s)$ and ${\mathcal
O}(\alpha_s^2)$ corrections to the $\chi_{c2}+\gamma$ production cross section
are sizable, which even reduces the Born order cross section by one order of
magnitude after including the NNLO perturbative correction. Taking the values
of the long-distance NRQCD matrix elements from nonrelativistic potential
model, our prediction to $\chi_{c1}+\gamma$ production rate is consistent with
the recent {\tt Belle} measurement. The NNLO predictions for the
$\chi_{c0,2}+\gamma$ production rates are much smaller than that for
$\chi_{c1}+\gamma$, which seems to naturally explain why the $e^+e^-\to
\chi_{c0,2}+\gamma$ channels have escaped experimental detection to date.
