  Inspired by precision tests of the Standard Model in future lepton colliders,
the numerical analysis of the following scattering processes, $e^+e^-
\rightarrow Z h^0\gamma$ and $e^+e^- \rightarrow Z H^0 \gamma$, are carried at
the tree level including all possible diagrams in the two-Higgs-doublet model
(2HDM). This model has many free parameters, but the parameters which take part
in the scattering amplitudes of these two processes are primarily the mixing
angle parameter, $s_{\beta-\alpha}$, and the masses of the neutral Higgs
bosons, ($h^0, H^0$). Therefore, measuring the production rates of $Z
h^0\gamma$ and $Z H^0 \gamma$ final states open another test for the scalar
sectors of the 2HDM. The numerical analysis is performed under the current
experimental constraints. The production rates and the asymmetry in the
forward-backward direction are presented as a function of the center-of-mass
energy covering the future lepton colliders. The unpolarized cross-section gets
up to 6.19 (4.86) fb at $\sqrt{s} = 350$ (500) GeV and 0.164 (0.157) fb at
$\sqrt{s} = 350$ (500) GeV for $e^+e^- \rightarrow Z h^0\gamma$ and $e^+e^-
\rightarrow Z H^0 \gamma$, respectively. The polarization of the incoming
$e^+e^-$ beams are studied for various configurations, and it enhances the
cross-section by a factor of 1.78 in both processes for
$P_{e^+,e^-}=(+0.6,-0.8)$.
