Fragmentation is the dominant production mechanism for heavy-quark-antiquark bound states with large transverse momentum. We analytically calculate the initial scale fragmentation functions (FFs) for a gluon to split into S-wave heavy quarkonium state using the perturbative QCD. Our analytical expression of FF depends on the transverse momentum k T of the gluon, and contains most of the kinematical and dynamical properties of the process. The analyses of this paper differ in that we present an analytical form of the transverse momentum dependent FF, using a different approach (Suzuki's model) in comparison with the numerical results presented in other references where the Braaten's model have been used. These k T dependent FFs are necessary to calculate the differential cross sections dσ /dk T , for which there are experimental data. We also incorporate, for the first time, hadron mass effects in our calculations to improve the FF at the small z regions, where z is the fragmentation parameter. These effects modify the relations between partonic and hadronic variables and reduce the available phase space and are responsible for the low-z threshold.