We establish the factorization formula for scalar Glueball production through radiative decays of vector states of heavy quarkonia, e.g. J/ψ, ψ(2S) and Υ(nS), where the Glueball mass is much less than the parent heavy quarkonium mass. The factorization is demonstrated explicitly at one-loop level through the next-to-leading order (NLO) corrections to the hard kernel, the non-relativistic QCD (NRQCD) long-distance matrix elements (LDMEs) of the heavy quarkonium, and the light-cone distribution amplitude (LCDA) of scalar Glueball. The factorization provides a comprehensive theoretical approach to investigate Glueball production in the radiative decays of vector states of heavy quarkonia and determine the physic nature of Glueball. We discuss the scale evolution equation of LCDA for scalar Glueball. In the end, we extract the value of the decay constant of Scalar Glueball from Lattice QCD calculation and analyze the mixing effect among f 0 (1370), f 0 (1500) and f 0 (1710).