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}$.