  Extensive numerical experiments on the long-term dynamics of planetesimals
near the orbits of planets around single stars with debris disks have been
carried out. The radial sizes of planetesimal clusters and the planetary
chaotic zone as a function of mass parameter $\mu$ (planet-star mass ratio)
have been determined numerically with a high accuracy separately for the outer
and inner parts of the chaotic zone. The results obtained have been analyzed
and interpreted in light of existing analytical theories (based on the
planet-planetesimal mean motion resonance overlap criterion) and in comparison
with previous numerical experiment approaches to the problem. We show and
explain how the stepwise dependence of the chaotic zone sizes on $\mu$ is
determined by the marginal resonances.
