  We calculate in the non-relativistic QCD (NRQCD) factorization framework all
leading relativistic corrections to the exclusive production of
$\chi_{cJ}+\gamma$ in $e^+ e^-$ annihilation. In particular, we compute for the
first time contributions induced by octet operators with a chromoelectric
field. The matching coefficients multiplying production long distance matrix
elements (LDMEs) are determined through perturbative matching between QCD and
NRQCD at the amplitude level. Technical challenges encountered in the
non-relativistic expansion of the QCD amplitudes are discussed in detail. The
main source of uncertainty comes from the not so well known LDMEs. Accounting
for it, we provide the following estimates for the production cross sections at
$\sqrt{s} = 10.6\textrm{ GeV}$: $\sigma (e^+ e^- \to \chi_{ c0} + \gamma) =
(1.3 \pm 0.4) \textrm{ fb}$, $\sigma (e^+ e^- \to \chi_{ c1} + \gamma) = (15.4
\pm 6.7) \textrm{ fb}$, and $\sigma (e^+ e^- \to \chi_{ c2} + \gamma) = (4.7
\pm 2.6) \textrm{ fb}$.
