Abstract: 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 μ (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 μ is determined by the marginal resonances.