Systematic uncertainties in the cross section measurement mainly come
from the luminosity measurement (\(\delta _{\cal L}\) ), branching fraction of \(h_c\rightarrow \gamma \eta _c\) , branching fractions of \(\eta _c\rightarrow X_i\) ,
detection efficiencies (\(\delta _{\epsilon _i\cdot {\cal B}(\eta _c\rightarrow X_i)}\) ),
radiative correction factors (\(\delta _{\rm ISR}\) ), vacuum polarization
factors (\(\delta _{\rm Vac}\) ) [1]}, and fits to
the mass spectrum. The integrated luminosity at each CME points is
measured using large-angle Bhabha events and has an estimated
uncertainty of 1.0%. The \(h_c\rightarrow \gamma \eta _c\) 
and \(\eta _c\rightarrow X_i\)  branching fractions are taken from Refs. [2]}, [3]}, and the
uncertainties in the radiative correction are the same as those used in the analysis of
\(e^+e^-\rightarrow \pi ^+\pi ^-h_c\)  [4]}. The uncertainties in the vacuum
polarization factor are 0.5% [1]}. The detection efficiency
uncertainty estimates are done with the same way as described in
Refs. [3]}, [7]}.
The uncertainty due to the \(\eta _c\)  mass (\(\delta _{\eta _c-\rm mass}\) ) is
estimated by changing its mass by \(\pm \) 1\(\sigma \)  of its world average value [8]};
the uncertainties due to the background shapes (\(\delta _{\rm bkg}\) ) are estimated
by changing the background function from a first-order to a second-order polynomial;
the uncertainty from the mass resolution (\(\delta _{\rm res}\) ) is
estimated by varying the mass resolution difference between data and MC simulation by
one standard deviation; the uncertainty from fit range (\(\delta _{\rm fit}\) ) is
estimated by extending the fit range; the uncertainty from the
\(\pi ^{0}\pi ^{0}h_{c}\)  substructure (\(\delta _{\rm sub}\) ) is estimated
by considering the efficiency with and without the inclusion of a
\(Z_{c}(4020)^{0}\) . The contribution from each source of systematic
error are listed in Table REF .
