We investigate the processes e+e− → ηc + γ and e+e− → χcJ + γ at B factories within the NRQCD factorization framework, computing the corresponding helicity amplitudes through \( \mathcal{O}\left({\alpha}_s^2\right) \) . The short-distance coefficients are obtained as series expansions in \( r=\frac{4{m}_c^2}{s} \) around r = 0, 1/3, 2/3, 1, using the method of differential equations. By combining the expansions from all four points, we construct composite asymptotic expressions that reproduce the exact results accurately over the full range 0 ≤ r ≤ 1, with relative errors below 0.1% over most of the domain and remaining under 1% elsewhere. Analytic expressions for the leading and next-to-leading logarithmic terms are extracted in the limit r → 0. Using these results, we compute the unpolarized cross sections and observe that the perturbative corrections are small for χc0 + γ, moderate for χc1 + γ, and substantial for ηc + γ and χc2 + γ. Theoretical prediction for χc1 + γ is consistent with the Belle measurement within 2σ, showing good agreement between theory and experiment. We also predict the angular distribution parameters \( {\alpha}_{\theta}^H \) , which are insensitive to NRQCD matrix elements and exhibit small theoretical uncertainties. These parameters further display good stability across different perturbative orders. With the high luminosity anticipated at Belle 2, future experimental measurements will thus provide a clear test of NRQCD factorization.
