Within the framework of nonrelativistic quantum chromodynamics (NRQCD) factorization, we compute the \( \mathcal{O}\left({v}^4\right) \) relativistic corrections to the fragmentation of a heavy quark into the color-singlet \( {}^1{S}_0^{\left[1\right]} \) and \( {}^3{S}_1^{\left[1\right]} \) quarkonium states. Using the Collins-Soper definition of the fragmentation function, we reproduce the known \( \mathcal{O}\left({v}^2\right) \) results. We find that the \( \mathcal{O}\left({v}^4\right) \) correction gives a positive contribution relative to the leading order result over a wide range of the light-cone momentum fraction z, while its magnitude remains much smaller than that of the \( \mathcal{O}\left({v}^2\right) \) correction. This behavior indicates a good convergence of the NRQCD relativistic expansion in this process. We further extend the calculation to the fragmentation functions in the unequal-mass case at \( \mathcal{O}\left({v}^4\right) \) and obtain the corresponding analytical expressions.
