We calculate the next-to-next-to-leading-order (NNLO) perturbative corrections to P -wave quarkonia annihilation decay to two photons, in the framework of nonrelativistic QCD (NRQCD) factorization. The order-α 2 s short-distance coefficients associated with each helicity amplitude are presented in a semi-analytic form, including the "light-by-light" contributions. With substantial NNLO corrections, we find disquieting discrepancy when confronting our state-of-the-art predictions with the latest BESIII measurements, especially fail to account for the measured χc2 → γγ width. Incorporating the effects of spin-dependent forces would even exacerbate the situation, since it lifts the degeneracy between the nonperturbative NRQCD matrix elements of χc0 and χc2 toward the wrong direction. We also present the order-α 2 s predictions to χ b0,2 → γγ, which await the future experimental test. Charmonium decay has historically played an important role in establishing the asymptotic freedom of QCD, and served as a clean platform to probe the interplay between pertubative and nonperturbative dynamics [1,2]. Among them, the electromagnetic decay χ c0,2 → γγ provide a particularly interesting, and, rich testing ground of QCD [3,4]. In the past decades, these decay channels have been extensively studied from various theoretical angles, such as nonrelativistic potential model [5,6], relativistic quark model [7][8][9], Bethe-Salpeter approach [10], nonrelativistic QCD (NRQCD) factorization [11,12], as well as lattice QCD [13]. On the experimental side, they were previously measured by . BESIII experiment [15] has recently reported their high precision results, Γ γγ (χ c0 ) = (2.33 ± 0.20 ± 0.13 ± 0.17) keV, (1a) Γ γγ (χ c2 ) = (0.63 ± 0.04 ± 0.04 ± 0.04) keV. (1b)