  The local number variance associated with a spherical sampling window of
radius $R$ enables a classification of many-particle systems in $d$-dimensional
Euclidean space according to the degree to which large-scale density
fluctuations are suppressed, resulting in a demarcation between hyperuniform
and nonhyperuniform phyla. To better characterize density fluctuations, we
carry out an extensive study of higher-order moments, including the skewness
$\gamma_1(R)$, excess kurtosis $\gamma_2(R)$ and the corresponding probability
distribution function $P[N(R)]$ of a large family of models across the first
three space dimensions, including both hyperuniform and nonhyperuniform models.
Specifically, we derive explicit integral expressions for $\gamma_1(R)$ and
$\gamma_2(R)$ involving up to three- and four-body correlation functions,
respectively. We also derive rigorous bounds on $\gamma_1(R)$, $\gamma_2(R)$
and $P[N(R)]$. High-quality simulation data for these quantities are generated
for each model. We also ascertain the proximity of $P[N(R)]$ to the normal
distribution via a novel Gaussian distance metric $l_2(R)$. Among all models,
the convergence to a central limit theorem (CLT) is generally fastest for the
disordered hyperuniform processes. The convergence to a CLT is slower for
standard nonhyperuniform models, and slowest for the antihyperuniform model
studied here. We prove that one-dimensional hyperuniform systems of class I or
any $d$-dimensional lattice cannot obey a CLT. Remarkably, we discovered that
the gamma distribution provides a good approximation to $P[N(R)]$ for all
models that obey a CLT, enabling us to estimate the large-$R$ scalings of
$\gamma_1(R)$, $\gamma_2(R)$ and $l_2(R)$. For any $d$-dimensional model that
"decorrelates" or "correlates" with $d$, we elucidate why $P[N(R)]$
increasingly moves toward or away from Gaussian-like behavior, respectively.

Statistical Mechanics;Soft Condensed Matter