This paper presents the first particle-mesh N-body simulations of the gravitating vacuum model, in which quantum vacuum energy gravitates locally through the ansatz ρ_vac(ρ_m) = Λ₀ − α·ρ_m, modifying the Poisson equation to ∇²Φ = 4πG(1+2α)ρ_m in overdense regions while keeping the background Friedmann equation identical to ΛCDM. The code is verified against ΛCDM at α = 0, and simulations are run for five values of α at 32³ and 64³ resolution with convergence analysis. THE CENTRAL DISCOVERY: NONLINEAR SCREENING The simulations revealed an effect that was not predicted by the linear perturbation theory of Paper #8: nonlinear self-screening of the gravitational modification. Linear theory predicted that α = +0.03 would enhance σ₈ by +35% relative to ΛCDM. The N-body simulation measured only +9.3% — a factor of ~3× suppression. For α = −0.003, linear theory predicted −2.9%; the simulation gave −0.9%. This is not a numerical artifact or a parameter adjustment. It is a direct, unavoidable consequence of the model's own physics. The mechanism is simple and robust: the modification G_eff = G(1+2α) applies only in overdense regions (δ > 0). In the early universe, the density field is nearly homogeneous — almost every point is near δ ≈ 0, and the modification acts on all matter, just as linear theory assumes. But as structure grows, overdense regions collapse into filaments and halos that occupy an ever-shrinking fraction of the total volume, while voids — where gravity is standard, unmodified — expand to fill the majority of space. By z = 0, the gravitational modification is effectively locked inside dense structures that occupy perhaps 20–30% of the volume. The volume-averaged σ₈ statistic, which integrates over 8 h⁻¹ Mpc spheres containing a mixture of voids and halos, dilutes the modification by the large void fraction. The model screens itself. Nobody designed this into the ansatz. It emerged from the simulation. When a model makes an unexpected, nontrivial prediction that the author did not anticipate — that is a good sign. It means the model has internal structure beyond what was put in by hand. RESOLUTION OF THE JWST–S₈ DUALITY Paper #8 identified a fundamental problem: the sign of α that resolves the S₈ tension (α 0 (constant, unchanging). At high redshift (z > 5), the universe is nearly homogeneous. The modification (1+2α) applies to essentially all matter. Growth is enhanced close to the full linear-theory prediction: D/D_ΛCDM ≈ 1 + 3α at z = 10. This is the regime relevant for JWST: massive galaxies form more efficiently because gravity is stronger everywhere. At low redshift (z 0.2 h/Mpc (nonlinear regime). This specific shape — no turnover at high k — distinguishes the model from f(R) gravity (chameleon screening) and DGP (Vainshtein screening) and is testable with Euclid at sub-percent precision. — Convergence study at 32³ and 64³ demonstrating that the σ₈ ratio (the physically meaningful quantity) converges at the ~1.5% level between resolutions, even though absolute σ₈ depends strongly on resolution. — Press–Schechter halo mass function ratios showing mass-dependent enhancement/suppression of cluster abundance, testable with eROSITA, Euclid, and the Vera C. Rubin Observatory. — f·σ₈(z) comparison with BOSS, WiggleZ, 6dFGS, and VIPERS data (best fit α = −0.003, Δχ² = −0.16 vs ΛCDM, conditional on fixed Planck cosmology). — Phase transition parameter space (z_c, α_h) mapping — retained as an exploratory scenario, but with reduced motivation given the screening result. — Density field slices from actual simulations showing the cosmic web morphology under modified gravity. — Complete, runnable simulation code provided (Python/NumPy/SciPy), ready to be scaled to production resolution (512³–1024³) on HPC facilities. This paper directly responds to the referee's challenge from Paper #8: "Without N-body simulations or at minimum a halo model treatment, the claim that σ₈ is reduced cannot be evaluated." The simulations have been run. The result was unexpected — and it strengthens the model. Keywords: N-body simulations, vacuum energy, structure formation, S8 tension, nonlinear screening, particle-mesh method, dark energy, growth factor, halo mass function, modified gravity, power spectrum, JWST, cosmic voids Related identifiers: Paper #8: "Structure Growth in the Gravitating Vacuum Model" (predecessor) Monograph: "What If the Vacuum Gravitates Locally?" (full program) License: Creative Commons Attribution 4.0 International (CC BY 4.0) Paper #9 in the research program "What If the Vacuum Gravitates Locally?"