We present the calculation at next-to-leading order (NLO) in α s of the fragmentation function of a gluon into heavy quarkonium in the color-octet spin-singlet S-wave channel. To calculate the real NLO corrections, we adapt a subtraction scheme introduced by Frixione, Kunszt, and Signer.Ultraviolet and infrared divergences in the real NLO corrections are calculated analytically by evaluating the phase-space integrals of the subtraction terms using dimensional regularization. The subtracted phase-space integrals are then evaluated in 4 space-time dimensions. The divergences in the virtual NLO corrections are also calculated analytically. After renormalization, all the divergences cancel. The NLO corrections significantly increase the fragmentation probability for a gluon into the spin-singlet quarkonium states η c and η b . The sum extends over the types of partons i (gluons, quarks, and antiquarks). The symbol "⊗" in Eq. (1) represents a convolution integral over the longitudinal momentum fraction z of the hadron H relative to the parton. The pQCD cross section dσ for producing the parton i can be expanded in powers of α s (p T ), and it includes convolutions with parton distributions if the colliding particles are hadrons. The nonperturbative factor D i→H (z) is a fragmentation function that gives the probability distribution for z. We refer to Eq. (1) as the leading-power (LP) factorization formula. It was derived by Collins and Soper in 1981 for the case of a light hadron H at a transverse momentum satisfying p T Λ QCD [1].The LP factorization formula in Eq. (1) applies equally well to heavy quarkonium at a transverse momentum satisfying p T m, where m is the mass of the heavy quark.A proof of this factorization theorem that deals with issues specific to heavy quarkonium production was first sketched by Nayak, Qiu, and Sterman in 2005 [2]. The LP factorization formula gives the leading power in the expansion in powers of m/p T . In the case of cross sections summed over quarkonium spins, the corrections are suppressed by m 2 /p 2 T . The LP factorization formula has limited predictive power, because the nonperturbative factors D i→H (z) are functions of z that must be determined from experiment.In 1994, Bodwin, Braaten, and Lepage proposed the NRQCD factorization formula for cross sections for heavy quarkonium production [3]. It uses an effective field theory called nonrelativistic QCD to separate momentum scales of order m and larger from momentum scales of order mv and smaller, where v is the typical relative velocity of the Q orQ in the quarkonium. The theoretical status of the NRQCD factorization conjecture is discussed in Ref. [4]. The predictive power of the LP factorization formula for heavy quarkonium in Eq. (1) can be increased by applying the NRQCD factorization formula to the fragmentation