  We discuss a new set of $\sim$ 500 numerical n-body calculations designed to
constrain the masses and bulk densities of Styx, Nix, Kerberos, and Hydra.
Comparisons of different techniques for deriving the semimajor axis and
eccentricity of the four satellites favor methods relying on the theory of Lee
& Peale (2006), where satellite orbits are derived in the context of the
restricted three body problem (Pluto, Charon, and one massless satellite). In
each simulation, we adopt the nominal satellite masses derived in Kenyon &
Bromley (2019a), multiply the mass of at least one satellite by a numerical
factor $f \ge 1$, and establish whether the system ejects at least one
satellite on a time scale $\le$ 4.5 Gyr. When the total system mass is large
($f \gg 1$), ejections of Kerberos are more common. Systems with lower
satellite masses ($ f \approx$ 1) usually eject Styx. In these calculations,
Styx often `signals' an ejection by moving to higher orbital inclination long
before ejection; Kerberos rarely signals in a useful way. The n-body results
suggest that Styx and Kerberos are more likely to have bulk densities
comparable with water ice, $\rho_{SK} \lesssim$ 2 g cm$^{-3}$, than with rock.
A strong upper limit on the total system mass, $M_{SNKH} \lesssim 9.5 \times
10^{19}$ g, also places robust constraints on the average bulk density of the
four satellites, $\rho_{SNKH} \lesssim$ 1.4 g cm$^{-3}$. These limits support
models where the satellites grow out of icy material ejected during a major
impact on Pluto or Charon.
