  The $e^+e^- \to HZ$ process with the subsequent decay of the Higgs boson $H
\to Z\gamma$ is studied, where both $Z$ bosons are reconstructed in the final
states with two jets. The analysis is performed using Monte Carlo data samples
obtained with detailed ILD (The International Large Detector) detector
simulation assuming an integrated luminosity of 2 ab$^{-1}$, beam polarizations
of ${\cal{P}}_{e^-e^+} = (-0.8, +0.3)$, and center-of-mass energy of $\sqrt{s}$
= 250 GeV. The analysis is repeated for the case of two 0.9 ab$^{-1}$ data
samples with polarizations ${\cal{P}}_{e^-e^+} = (\mp0.8, \pm0.3)$.
Contributions of the potential background processes are studied using all
available ILD MC event samples. The largest background comes from the $e^+e^-
\to W^+W^-$ process supplemented by an energetic photon produced by initial
state radiation. To suppress this background we require that at least one of
the $Z$ bosons decays to $b$-jets. To reduce the jet reconstruction
uncertainties the $M_{\Delta} = M(jj\gamma) - M(jj) + M(Z_{\rm nom})$ variable
is used, where $M(Z_{\rm nom})$ = 91.2 GeV. The $M_{\Delta}$ distributions are
obtained for the studied signal and backgrounds to estimate the expected
accuracy of the ${\cal B}r(H \to Z\gamma)$ measurement. The accuracy is 22$\%$
for the option of the single polarization sample described above and
deteriorate to 24$\%$ in case of the sample with two polarizations. The
proposed method can be applied at any future $e^+e^-$ collider.
