  In this research, we conduct a global QCD analysis of fragmentation functions
(FFs) for neutral pions ($\pi^0$), neutral kaons ($K_S^0$), and eta mesons
($\eta$), utilizing world data of single inclusive hadron production in
$e^+e^-$ annihilation involving the most recent BESIII data with low collision
energy, to test the operational region of QCD collinear factorization for
single inclusive hadron production. We found that the QCD-based analysis at
next-to-next-to leading order in perturbative QCD with parameterized
higher-twist effects can explain both existing high-energy world data and the
BESIII new measurements, while the latter cannot be explained with existing FFs
extracted with high-energy data. To investigate the higher-twist contributions
to this discrepancy, a direct functional approach is employed, providing
testing framework for characterizing the experimental results over a wide range
of energy scales, from low to high, thus extending the classical theoretical
models to the BESIII domain.
