The next-to-next-to-leading order (NNLO) pQCD prediction for the $\gamma\gamma^* \to \eta_c$ form factor was evaluated using nonrelativistic QCD (NRQCD) in 2015. A strong discrepancy between the NRQCD prediction and the BaBar measurements was observed. Until now there has been no solution for this puzzle. In this paper, we present a NNLO analysis by using the Principle of Maximum Conformality (PMC). By carefully dealing with the light-by-light diagrams at the NNLO level, the resulting high precision PMC prediction agrees with the BaBar measurements within errors, and the conventional renormalization scale uncertainty is eliminated. The PMC is consistent with all of the requirements of the renormalization group, including scheme independence. The application of the PMC thus provides a rigorous solution for the $\gamma\gamma^* \to \eta_c$ form factor puzzle, emphasizing the importance of a correct renormalization scale-setting. The results also support the applicability of NRQCD to hard exclusive processes involving charmonium. Show Partial Abstract