Based on the NRQCD factorization formalism, we calculate the next-to-nextto-leading order QCD corrections to the heavy quarkonium η c (η b ) production associated with a photon at electron-positron colliders. By matching the amplitudes calculated in full QCD theory to a series of operators in NRQCD, the short-distance coefficients up to NNLO QCD radiative corrections are determined. It turns out that the full set of master integrals that we obtained could be analytically expressed in terms of Goncharov Polylogarithms, integrals over polylogarithms and complete elliptic integrals, which mostly do not exist in the literature and could be employed in the analyses of other physical processes. In phenomenology, numerical calculations of NNLO K-factors and cross sections of e + e − → γ + η c (η b ) processes in BESIII and B-factory experiments are performed, which may stand as a test of the NRQCD higher order calculation while confronting to the data.