An unbinned maximum likelihood fit is applied to the \(M^{\rm recoil}_{\pi ^0}|_{\rm max}\) 
distribution summed over all 16 \(\eta _c\)  decay modes. The data at
\(\sqrt{s}=4.23\) , 4.26, and 4.36 GeV are fitted simultaneously with
the same signal function with common mass and width.
The signal shape is parametrized with a constant-width
relativistic Breit-Wigner function convolved with a Gaussian-distributed
mass resolution, where
the mass resolution is determined from a fit to a MC sample with the
width set to zero. Because of the limited statistics of the \(Z_c(4020)^0\) 
signal,
its width is fixed to that of its charged partner, (\(7.9\pm 2.6\) ) MeV [1]}.
Assuming the spin and parity of the \(Z_c(4020)^0\)  are
\(1^+\) , a phase space factor \(pq^3\)  is included in the partial width,
where \(p\)  is the \(Z_c(4020)^0\)  momentum in the \(e^+e^-\)  rest frame and \(q\)  is
the \(h_c\)  momentum in the \(Z_c(4020)^0\)  rest frame.
