Local Number Fluctuations in Hyperuniform and Nonhyperuniform Systems: Higher-Order Moments and Distribution Functions

The local number variance σ(R) associated with a spherical sampling window of radius R enables a classification of many-particle systems in d-dimensional Euclidean space R according to the degree to which large-scale density fluctuations are suppressed, resulting in a demarcation between hyperuniform and nonhyperuniform phyla. To more completely characterize density fluctuations, we carry out an extensive study of higher-order moments or cumulants, including the skewness γ1(R), excess kurtosis γ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 systems with varying degrees of shortand long-range order. To carry out this comprehensive program, we derive new theoretical results that apply to general point processes and conduct high-precision numerical studies. Specifically, we derive explicit closed-form integral expressions for γ1(R) and γ2(R) that encode structural information up to three-body and four-body correlation functions, respectively. We also derive rigorous bounds on γ1(R), γ2(R) and P [N(R)] for general point processes and corresponding exact results for general packings of identical spheres. High-quality simulation data for γ1(R), γ2(R) and P [N(R)] are generated for each model. We also ascertain the proximity of P [N(R)] to the normal distribution via a novel Gaussian “distance” metric l2(R). Among all models, the convergence to a central limit theorem (CLT) is generally fastest for the disordered hyperuniform processes such that γ1(R) ∼ l2(R) ∼ R−(d+1)/2 and γ2(R) ∼ R−(d+1) for large R. 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 across all dimensions for intermediate to large values of R, enabling us to estimate the large-R scalings of γ1(R), γ2(R) and l2(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. Our work elucidates the fundamental importance of higher-order structural information to fully characterize density fluctuations in many-body systems across length scales and dimensions, and thus has broad implications for condensed matter physics, engineering, mathematics and biology.

Consider a statistically homogeneous (translationally invariant) point process in d-dimensional Euclidean space R d at number density ρ and sampling for the number of points N (R) within a d-dimensional spherical window of radius R (see Fig. 1) and volume The local number variance σ 2 (R) ≡ N 2 (R) − N (R) 2 is a useful measure of number fluctuations, where the first moment N (R) = ρv 1 (R) is the average number of points within a d-dimensional spherical (sampling) window of radius R and angular brackets denote an ensemble average. The local number variance is exactly determined by pair statistics, and can be given either in terms of the pair correlation function g 2 (r) in direct space or the structure factor S(k) in reciprocal space [13]: = ρv 1 (R) 1 (2π) d where h(r) ≡ g 2 (r) − 1 is the total correlation function, α 2 (r; R) is the intersection volume of two spherical windows of radius R, scaled by v 1 (R), whose centers are separated by the distance r, andα 2 (k; R) is its Fourier transform. The large-scale behavior of the number variance σ 2 (R) is central to the hyperuniformity concept, which is attracting attention across many fields [13,16,[19][20][21][22][23]. Specifically, a hyperuniform point process is one in which σ 2 (R) grows slower than the window volume, i.e., R d , for large R and hence is characterized by large-scale density fluctuations that are anomalously suppressed compared to those of typical disordered systems. The hyperuniformity concept generalizes the traditional notion of long-range order of crystals and quasicrystals to also encompass certain exotic disordered states of matter [13,16]. Disordered hyperuniform systems are diametrically opposite to systems at thermal critical points in which the local variance diverges faster than R d in the limit R → ∞. Any system with such divergent behavior in the local variance has been called anti-hyperuniform [16] (see also Sec. II for additional details).
While the local number variance contains useful information, one would like to more completely characterize the fluctuations by ascertaining higher-order moments, i.e., N k (R) , where k ≥ 3, as well as the corresponding discrete probability distribution P [N (R)] associated with finding exactly N (R) particles within a d-dimensional spherical window of radius R. The mth moment of the distribution is given by Due to the fact that the random variable N (R) is discrete and cannot take on negative values, the probability distribution P [N (R)], for finite R, can never exactly attain the normal distribution, which is given by For example, for a statistically homogeneous Poisson point process in R d at number density ρ, which deviates significantly from the normal distribution for sufficiently small R. This is one of the rare cases in which a closed-form analytic formula for P [N (R)] is known across dimensions for nontrivial point processes.
It is only when R tends to infinity that the Poisson distribution becomes a normal distribution, i.e., it follows a central limit theorem (CLT); see Ref. [24] and references therein. The reader is referred to Refs. [25][26][27][28] for proofs of CLT for other point processes. In this paper, we investigate the skewness γ 1 (R) (related to the first three moments), excess kurtosis γ 2 (R) (related to the first four moments), and the number distribution function P [N (R)] for general homogeneous point processes in R d as well as a wide class of models across dimensions. Specifically, we derive explicit closedform integral expressions for γ 1 (R) in terms of the number density ρ, the pair correlation function g 2 and the three-body correlation function g 3 (defined in Sec. II). Similarly, we derive corresponding formulas for the excess kurtosis γ 2 (R), which now depends additionally on the four-body correlation function g 4 . These integral relations also involve geometrical information about the spherical windows via the intersection volumes of up to three and four spheres in the cases of the skewness and excess kurtosis, respectively. Thus, the skewness and excess kurtosis encode up to three-body and four-body information about spatial correlations and window geometries, respectively. We also derive some exact elementary results for the skewness, excess kurtosis and number distribution that apply to general packings of identical spheres.
Via high-precision computer simulation studies, we accurately determine the number variance, skewness, excess kurtosis and the number distribution for up to eight different models of statistically homogeneous point processes across the first three space dimensions and wide range of window radii R. These models include both nonhyperuniform and hyperuniform systems with varying degrees of short-and long-range order. Among the nonhyperuniform point processes, we characterize fluctuations of Poisson, random sequential addition (RSA) packings, equilibrium hard spheres, Poisson cluster and anti-hyperuniform point processes. Among the hyperuni-form point processes, we characterize hypercubic lattices, randomly perturbed lattices, and stealthy disordered hyperuniform systems. We show that our simulation results for γ 1 (R), γ 2 (R) and P [N (R)] for all models are in excellent agreement with the aforementioned rigorous bounds and exact results for the applicable ranges of R. For all disordered hyperuniform models, our explicit general formulas of these quantities in terms of n-body information enable us to infer the existence of "hidden" order that manifests itself for the first time at the three-body level or higher.
For each model considered in this paper, we are interested in ascertaining how large R must be such that P [N (R)] is well-approximated by the Gaussian (normal) distribution. We have found that such "distance" metrics proposed previously are not adequate for assessing the diverse set of models that we consider here across dimensions. We quantify this proximity to the normal distribution for any model by introducing a certain Gaussian "distance" metric l 2 (R), defined in Sec. VII. This distance metric enables us to accurately determine when the distribution function P [N (R)] for a particular model is tending to a CLT. Because the distributions for all models (except the lattices) across dimensions are unimodal, the tendency to a CLT corresponds to the skewness and excess kurtosis simultaneously tending to zero. We have found that almost all of the considered models across dimensions obey a CLT. The convergence to a CLT is slowest for the anti-hyperuniform point process, followed by the Poisson cluster process, and the Poisson process. The nonhyperuniform RSA and equilibrium packings tend to a CLT at the same rate as a Poisson point process but with smaller coefficients of proportionality. Among all models, the convergence to a CLT is generally fastest for the disordered hyperuniform processes. The only models considered that do not achieve a CLT are the hypercubic lattices for any d and 1D hyperuniform systems of class I. The reader is referred to Sec. VII for details.
We have examined a variety of well-known closedform probability distributions P [N (R)] to ascertain those which best approximate the actual distributions for finite R for all of our models. Interestingly, the gamma distribution provides a good approximation to the number distribution P [N (R)] for all models that obey a CLT across all dimensions for intermediate to large values of R (Sec. VII C). It is noteworthy that the approximation of P [N (R)] by a gamma distribution enables us to estimate the large-R scalings of γ 1 (R), γ 2 (R) and l 2 (R) for all models across dimensions that obey a CLT. Among all models, the convergence to a CLT is generally fastest for the disordered hyperuniform processes such that γ 1 (R) ∼ l 2 (R) ∼ R −(d+1)/2 and γ 2 (R) ∼ R −(d+1) for large R. For standard nonhyperuniform models, convergence to a CLT is slower such that γ 1 (R) ∼ l 2 (R) ∼ R −d/2 and γ 2 (R) ∼ R −d . Finally, convergence to a CLT is slowest for the antihyperuniform model such that γ 1 (R) ∼ l 2 (R) ∼ R −1/2 and γ 2 (R) ∼ R −1 . These predictions are corroborated by corresponding simulation re-sults.
An important fundamental question is what is the effect of increasing the space dimension on number fluctuations for any particular model? To answer this question, we recall the so-called decorrelation principle [29], which roughly states that for any disordered point process, unconstrained correlations that exist in low dimensions vanish as d tends to infinity, and all higher-order correlation functions g n for n ≥ 3 may be expressed in terms of the number density ρ and pair correlation function g 2 . The decorrelation principle was employed to justify the conjecture that the densest sphere packings in sufficiently high dimensions are disordered (as opposed to ordered in low dimensions) [29,30]. Importantly, decorrelation in pair statistics has been shown to manifest itself in low dimensions in the case of disordered sphere packings [31][32][33] as well as other disordered systems with strongly repulsively interacting particles [34,35]. Since the number distribution function P [N (R)] generally involves certain integrals over all of the n-body correlation functions, the decorrelation principle implies that P [N (R)] increasingly becomes Gaussian-like as the space dimension increases for any model that decorrelates with d. Similarly, for models the correlate with d, P [N (R)] increasingly deviates from the normal distribution as d increases. We confirm such behaviors for all models that obey a CLT (Sec. VII D).
In Sec. II, we provide basic definitions and necessary background material. In Sec. III, we derive explicit closed-form integral expressions for γ 1 (R) and γ 2 (R) as well as rigorous lower bounds on both of these quantities. We also obtain some general exact results for the first few cumulants and distribution functions for sphere packings, whether disordered or not. Section IV describes the large variety of nonhyperuniform and hyperuniform models in one, two and three dimensions that we study in this paper. In Sec. V, we discuss our proposed Gaussian distance metric. Section VI describes the simulation procedure that we employ to sample the first four cumulants and number distributions. In Sec. VII, we present our results. We make concluding remarks in Sec. VIII.

II. DEFINITIONS AND BACKGROUND
A stochastic point process in R d is defined as a mapping from a probability space to configurations of points r 1 , r 2 , r 3 . . . in d-dimensional Euclidean space R d ; see Ref. [36] for mathematical details. Let X denote the set of configurations such that each configuration x ∈ X is a subset of R d that satisfies two regularity conditions: (i) there are no multiple points (r i = r j if i = j) and (ii) each bounded subset of R d must contain only a finite number of points of x (i.e.,, x is "locally finite"). The point process is statistically is characterized by the generic nparticle probability density function ρ n (r n ), where r n is a shorthand notation for the position vectors of any n points, i.e., r n ≡ r 1 , r 2 , . . . , r n [5,37]. In words, the quantity ρ n (r n )dr n is proportional to the probability of finding any n particles with configuration r n in volume element dr n ≡ dr 1 dr 2 · · · dr n , i.e., it is the probability measure. For any subvolume Ω ∈ R d , the following normalization (average) condition involving the fluctuating number of particles within this subvolume, N Ω , immediately follows Note that this random setting is quite general; it incorporates cases in which the locations of the points are deterministically known, such as a lattice. For statistically homogeneous media, ρ n (r n ) is translationally invariant and hence depends only on the relative displacements, say with respect to r 1 : ρ n (r n ) = ρ n (r 12 , r 13 , . . . , r 1n ), where r ij = r j − r i . In particular, the one-particle function ρ 1 is just equal to the constant number density of particles ρ. For statistically homogeneous point patterns, it is convenient to define the so-called n-particle correlation function In systems without long-range order and in which the particles are mutually far from one another, ρ n (r n ) → ρ n and we have from (9) that g n (r n ) → 1. Thus, the deviation of g n from unity provides a measure of the degree of spatial correlation between the particles, with unity corresponding to no spatial correlation. The important two-particle quantity g 2 (r 12 ) is usually referred to as the pair correlation function. The total correlation function h(r 12 ) is defined as h(r 12 ) = g 2 (r 12 ) − 1, and thus is a function that is zero when there are no spatial correlations in the system. Observe that the structure factor S(k) is related to the Fourier transform of h(r), denoted byh(k), via the expression whereh A hyperuniform point process is one in which singlescattering events at infinite wavelength vanishes, i.e., lim |k|→0 S(k) = 0.
This implies that a hyperuniform system obeys the following sum rule: in direct space and hence h(r) must exhibit negative pair correlations, i.e., anticorrelations, for some values of r [16]. By contrast, an anti-hyperuniform point process is one in which S(k) tends to +∞ in the limit |k| → 0 [16]. A lattice Λ in R d is a subgroup consisting of the integer linear combinations of independent vectors that span R d and thus represents a special subset of point processes. In a lattice Λ, the space R d can be geometrically divided into identical regions F called fundamental cells, each of which contains just one point specified by the lattice vector [38,39]. In the physical sciences, a lattice is equivalent to a Bravais lattice. Unless otherwise stated, we will use the term lattice. A periodic point process is a more general notion than a lattice because it is is obtained by placing a fixed configuration of N points (where N ≥ 1), called the basis, within one fundamental cell of a lattice Λ, which is then periodically replicated. Thus, the point process is still periodic under translations by Λ, but the N points can occur anywhere in the chosen fundamental cell. Any lattice or periodic point configuration can be made statistically homogeneous by uniform translations of the pattern within the fundamental cell.
We call a packing in R d a collection of nonoverlapping particles [40]. The centroids of the particles constitute a special point process in which no two particles can closer than some minimal distance. In this paper, we consider packings of identical spheres of diameter D. The packing fraction φ = ρv 1 (D/2) is the fraction of space covered by the spheres, where v 1 (R) is the volume of a sphere of radius R given by (1). Any periodic point configuration with a finite basis can be regarded to be a packing since there is a minimal pair distance.

III. GENERAL LOCAL MOMENT FORMULAS AND PROBABILITY DISTRIBUTION FOR HOMOGENEOUS POINT PROCESSES
We consider the determination of local moment formulas using the formalism of Torquato and Stillinger [13] that was used to obtain formulas for the local number variance. Here we immediately begin with a d-dimensional spherical window of radius R in ddimensional Euclidean space R d with window indicator function where x is some arbitrary position vector in R d and x 0 is the position vector of the window center. The number of points N (R; x 0 ) within the window at position x 0 is given by which must be a finite number. We will subsequently use the fact that v int n (r n ; is the intersection volume of n spheres of radius R centered at positions r 1 , r 2 , · · · , r n .

A. Moments
The mth moment associated with the random variable N (R) is given by the following ensemble average: Here we have implicitly assumed homogeneity, which renders the mth moment independent of the position of the window x 0 . Under the ergodic assumption, the ensemble average indicated on the left-hand side of relation (18) is equivalent to averaging by uniformly window sampling a single realization over the infinite space.
Following the same procedure used in Ref. [13] to obtain an explicit formula for the second moment, we obtain from (18) that the mth moment N m (R) for a homogeneous process is given by integrals involving the finite set of correlation functions g 2 , g 3 , . . . g m weighted with the set of intersection volumes v int 2 , v int 3 , . . . , v int m . For example, expanding the sums in (18) into one-body, two-body, . . . and m-body terms in a manner analogous to the one given in Ref. [13], third and fourth moments are explicitly given by and where v int n is given by (17).
We are generally interested in the mth order cumulant C m (R), which is directly related to the mth central moment [N (R) − N (R) ] m . For example, the first several cumulants are given by and where we have used the following identities: We now derive a lower bound on C 3 (R) in terms of the first and second cumulants for general point processes. This easily follows from the fact that the product g 3 (r 3 )v int 3 (r 3 ; R) in the last integral of (22) is nonnegative for all positions and so dropping this integral yields the following lower bound: Similarly, dropping the positive integral in (23) involving the product g 4 (r 4 )v int 4 (r 4 ; R) gives the following lower bound on C 4 (R) in terms of the first, second and third cumulants: These bounds are relatively tight for sufficiently small values of R and can be exact (sharp) for such R for packings, as discussed in Sec. III C.
The third-and fourth-order cumulants are directly related to the skewness and excess kurtosis, respectively. The skewness is often defined as which, qualitatively speaking, is a measure of the asymmetry of the probability distribution. The excess kurtosis, which is a measure of the heaviness of the "tails" of the probability distribution, is defined to be For a general random variable, Pearson derived a lower bound on the excess kurtosis in terms of the skewness [41]: We also apply this lower bound to validate our numerical results for all of our model point processes. Such details are reported in the SM. Both γ 1 and γ 2 are identically zero for the normal distribution and hence such lower-order information can herald at what value of R a general point configuration can be approximated by a normal distribution. Note that the higher-order cumulants become increasingly more complicated but can be written as a determinant involving the moments N m (R) [42].
In the case of a homogeneous Poisson point process, g n = 1 for all n. Hence, from the expressions above, it immediately follows that C 2 (R) = C 3 (R) = C 4 (R) = ρv 1 (R), which are the expected well-known results for this point process. It follows that the corresponding skewness and excess kurtosis are exactly given by and respectively. Observe that both γ 1 (R) and γ 2 (R) tend to zero in the limit R → ∞, which is consistent with the fact that the Poisson point process obeys a CLT.

B. Probability Distribution Function
It is straightforward to show that for an arbitrarilyshaped window region Ω, the probability distribution function is given by [2,3] where (37) is exactly the same as the average given in (7) under the assumption of homogeneity. Thus, we see that the average A m (Ω) is the nontrivial and common contribution to P [N Ω ] for any specific value of N Ω . The only differences in P [N Ω ] for different values of N Ω are the combinatoric factors multiplying the coefficients A m (Ω) in the series (36). Because these coefficients are intrinsically positive, relation (36) is an alternating series.
When Ω is a spherical window of radius R, it simply follows that where For a fixed value of M ∈ N 0 , P [N (R) < M ] as a function of R is a complementary cumulative distribution function, which is associated with the "void" probability density function H V (R; M ) [9], where H V (R; M )dR is the probability that the distance to the M th nearest neighbor from an arbitrary point in space is between R and R + dR. For example, in the case M = 1, , which is associated with the void nearest-neighbor probability density function H V (R) [37,39,43]. In the case of a homogeneous Poisson point process, we can immediately recover the exact result (6) for the number distribution from relations (38) and (39) using the fact that g n = 1 for all n and the identity Exact results for P [N (R)] for non-Poissonian point processes are rare. One exception is the one-dimensional model of equilibrium hard rods for which P [N (R)] is known exactly [9]. In the case of disordered equilibrium hard-disk (d = 2) and hard-sphere (d = 3) packings, accurate but approximate expressions for the exclusion probability E V (R) ≡ P [N (R) = 0] are available [37,43]. Thus, in principle, one can extract from such formulas approximations for A m (R) to get corresponding approximations for P [N (R)] for any N (R) ≥ 1 for such packings. The difficulty in ascertaining P [N (R)] exactly for nontrivial models can be appreciated by appealing to the ghost random sequential addition (RSA) packing process [29], for which the g m are known exactly for any m. While the evaluation of the integral (39) for ghost RSA can be carried out exactly for very small m, its exact determination becomes impossible for general values of m. This points to the importance of devising accurate numerical methods to determine the number distributions of packings. For example, P [N (R)] has be determined in simulations for equilibrium hard spheres and MRJ sphere packings in Ref. [15].
It has already been established rigorously that truncations of the alternating series (38) for the special case P [N (R) = 0] = E V (R) at an even and odd number of terms yield successive upper and lower bounds on the probability distribution P [N (R) = 0], respectively. Such bounds are consequences of an inclusion-exclusion principle associated with this alternating series. Here we make the simple observation that the same inclusion-exclusion principle applies for any value of N (R). The first several of such bounds are given by where A N (R) is a shorthand for A N (R) (R). These bounds become increasingly sharper as more terms are included. Moreover, these bounds can be sharp (exact) for sufficiently small R for sphere packings, as discussed in Sec. III C. We utilize these bounds to validate our simulations results for all models across dimensions in the SM.

C. Elementary Results for Packings
Here we obtain some general exact results for the first few cumulants and distribution functions for packings of identical spheres of diameter D, whether disordered or not. Results that apply to lattice packings are also derived.
Because no two spheres can overlap when their centers are separated by a distance less than or equal to D, the cumulants can be written explicitly for d ≥ 1 and R ≤ D/2, since such a window region can accommodate at most a single sphere. For example, this means that integral involving g 2 in the last line of Eq. (21) is identi-cally zero, and hence for R ≤ D/2, which was noted by Torquato and Stillinger [13]. More generally, noting that N (R) m = ρv 1 (R) for all m because all terms in (19) and (20) involving either g 2 , g 3 and g 4 , we have from (22) and (23) for d ≥ 1 that for R ≤ D/2, This means that for R ≤ D/2, and As in the Poisson case, both the skewness γ 1 (R) and excess kurtosis γ 2 (R) diverge to +∞ in the limit R → 0. For 0 ≤ R ≤ D/2, while γ 1 (R) is generally a monotonically decreasing function of R that can be nonnegative, γ 2 (R) is generally a nonmonotonic function but also can be nonnegative. The corresponding results for the distribution function P [N (R)] for R ≤ D/2 follow immediately from (36) and (39), namely, When D/2 < R ≤ D/ √ 3 and d ≥ 2, the window can accommodate at most two but not three hard spheres. Therefore, following the same reasoning as above for R ≤ D/2, we find from (22) and (23) that for d ≥ 2 and and We see that both C 3 (R) and C 4 (R) are given purely in terms of the mean N (R) = ρv 1 (R) and number variance σ 2 (R) for such R. Observe also that formula (53) is identical the lower bound (29) for general point processes. It is seen that for R ≤ D/2, relations (46) and (47) are recovered, as expected. Finally, for d = 1, these formulas actually apply for R ≤ D.
When the window can accommodate at most three spheres, the fourth cumulant C 4 (R) for 0 ≤ R ≤ R * (d) is exactly given in terms of the first three cumulants, i.e., where R * (d) > D/ √ 3 is a threshold that depends on the space dimension d. For example, R * (1) = 3D/2, R * (2) = √ 2D/2, and R * (3) = 3/8D. Note that formula (55) is identical to the lower bound (30) for general point processes.
When the spherical window can accommodate at most two spheres, implying that R ≤ D/ √ 3 for d ≥ 2, three-body and higher-order terms in the series (39) for the probability distribution vanish identically, yielding the following exact result for P [N (R)]: We see that for such windows, the entire distribution function is completely determined by the first and second cumulants. A nontrivial upper bound on the variance σ 2 (R) of a packing follows immediately from (57) and the fact that P [N (R) = 1] must be nonnegative, i.e., For the same reasons, relation (58) yields the following general lower bound on the variance As before, these formulas actually apply for R ≤ D for d = 1.
Similarly, when the spherical window can accommodate at most three spheres, implying that 0 ≤ R ≤ R * (d), four-body and higher-order terms in the series (39) for P [N (R)] vanish identically, yielding the following exact result for P [N (R)]: Thus, for such R, the entire probability distribution function is completely determined by the first three cumulants.  65)] yield lower bounds on C 3 (R) in terms of the first two cumulants and so the maximum of these two lower bounds are to be chosen. In the SM, we demonstrate that the aforementioned exact results for C 3 (R), C 4 (R) and certain P [N (R)] are in excellent agreement with our corresponding simulations data for sphere packings examined in the paper across dimensions.
More generally, for any packing of identical spheres, a spherical region of radius R can accommodate a maximum number of spheres, denoted by N max (R). This maximal number can be determined from tabulations of the so-called densest local packings for a finite range of particle numbers in both two [44] and three [45] dimen-sions. Therefore, the number distribution P [N (R)] for a packing is generally far from a normal distribution for a finite-sized window, since it must have compact support such that it is zero for N (R) > N max (R), i.e., In such instances, the distribution function is determined by a finite set of moments, i.e., the first, second, . . ., N max th moments. Moreover, for any dense packing or point process in which the nearest neighbor from a particle is narrowly distributed (e.g., "strongly" stealthy systems described below and in Sec. IV B 3), P [N (R)] will be zero for N (R) < N min (R), where the cut-off value N min (R) grows with R. This situation prevents a strict CLT from applying for finite-sized windows.
Another important observation is that for point processes in which the "hole" radius R is bounded from above by R max , the probability of finding a spherical window with radius R > R max must be zero, i.e., P [N (R) = 0] = 0 for R > R max , which of course is non-Gaussian behavior. The cut-off value R max for a point process in R d is its covering radius [39]. Processes with this bounded-hole property include periodic packings with a finite basis [46], quasicrystals [47], as well as the saturated random sequential addition packing process (see Sec. IV A 3). Disordered stealthy point processes have bounded holes [46,48], as discussed in Sec. IV B 3.
We note that for the hypercubic lattice Z d scaled by D (see Sec. IV B 1 for precise definition), all of the relations derived above for the skewness, excess kurtosis and distribution function for the situation R ≤ D/ √ 3 actually apply as well for the larger range R ≤ D/ √ 2 when d ≥ 2, where D is the lattice spacing. In fact, in the case of the scaled integer lattice ρ −1 Z at number density ρ, we can obtain an exact formula for P [N (R)] by invoking the key idea of Ref. [13] to yield the exact result for the local number variance, namely, the number of points inside a window of radius R can only take two values, either N R or N R + 1, where N R ≡ ρ2R and x is the floor function of a real number x. The probability distribution for all N (R) is given by where {x} ≡ x − x is the fractional part of a positive number x. Thus, this skewed distribution is highly non-Gaussian with nonexistent left or right tails for almost all values of N (R), implying values of the skewness and excess kurtosis that are generally far from zero for almost R. From the distribution function (71) and relation (4), we can immediately obtain the first several cumulants: which are all periodic functions with period ρ2R. Relation (72) for the variance was given in Ref. [13]. The reader is referred to the top panel of Fig. 2, which shows plots of the variance, skewness and excess kurtosis as a function of R for Z. The highly discrete nature of the number distribution for the integer lattice extends to that for the hypercubic lattice Z d for d ≥ 2, as will see in Sec. VII.

IV. NONHYPERUNIFORM AND HYPERUNIFORM MODELS
We consider eight different models of statistically homogeneous point processes in two and three dimensions: five nonhyperuniform models, one of which is antihyperuniform (hyperplanes intersection process or HIP), and three hyperuniform models. Analogous models are also examined in one dimension, except for HIP, which is not defined in this dimension. The reader is referred to Figs. 3 and 4 for representative images of configurations for each of the models in two dimensions.
It is useful to recall scaling relations for hyperuniform and nonhyperuniform point processes. Consider any homogeneous point process in R d for which the structure factor has the following power-law behavior as the wavenumber tends to zero: This scaling implies that the total correlation function h(r) has the corresponding power-law behavior 1/|r| d+α for large |r| [16]. For hyperuniform systems, the exponent α is a positive constant, which implies that there are three different scaling regimes (classes) that describe the associated large-R of the number variance [13,16,49]: By contrast, for any nonhyperuniform system, it follows from the asymptotic analysis given in Ref. [16] that The scaling for the anti-hyperuniform instance can be obtained using an asymptotic analysis of either the directspace representation (2) or the Fourier-space representation (3) of the number variance, as derived in Ref.
[50]. The typical nonhyperuniform scaling in (77) results from the fact that S(0) is bounded and, indeed, the implied constant multiplying R d is proportional to S(0). Any nonhyperuniform point process for which S(0) > 1 has a large-R asymptotic number variance σ 2 (R) that is larger than that for a Poisson point process [S(0) = 1] with the same mean N (R) . We call such a nonhyperuniform point process super-Poissonian. Two examples of super-Poissonian point processes studied in this work are the Poisson cluster and HIP point processes described below.
A. Nonhyperuniform Processes

Poisson Point Process
A homogeneous Poisson point process in R d has a structure factor S(k) = 1 for all k and hence is nonhyperuniform. At unit mean density (ρ = 1) this process is generated within a hypercubic simulation box of fixed volume V under periodic boundary conditions by a twostep procedure. First, we choose a random number N from the Poisson distribution (6) with intensity or mean ρV = V and then place N points in the simulation box uniformly.

Equilibrium Packings
We also consider equilibrium packings of identical sphere (Gibbs hard-sphere processes) across the first three space dimensions. For d = 2 and d = 3, we examine disordered states that lie along the stable liquid branch [5,37] as well as disordered states in one dimension, all of which are nonhyperuniform. We generate such equilibrium packings using the well-established Metropolis numerical scheme [5,37]. All configurations that we generate are well away from jamming points and hence all are nonhyperuniform with bounded S(0) (see Table  I).

Random Sequential Addition Packings
The random sequential addition (RSA) process is a time-dependent (nonequilibrium) procedure that generates disordered sphere packings in R d [32,[51][52][53][54][55]. Starting with an empty but large volume in R d , the RSA process is produced by randomly, irreversibly, and sequentially placing nonoverlapping spheres into the volume. This procedure is repeated for ever increasing volumes and then an appropriate infinite-volume limit is obtained. In practice, hard spheres are randomly and sequentially placed into a large fundamental cell under periodic boundary conditions and subject to a nonoverlap constraint: If a new sphere does not overlap with any existing spheres, it will be added to the configuration; otherwise, the attempt is discarded. One can stop the addition process at any time t, obtaining RSA configurations with a range of packing fractions φ(t) up to the maximal saturation value φ s ≡ φ(∞), which imposes a bounded-hole property [46]. For identical spheres, which we consider here, φ s ≈ 0.74, 0.55, and 0.38 for d = 1, 2, and 3, respectively [33,51,[53][54][55]. The pair correlation function g 2 (r) is known exactly only in one dimension [56]. The structure factors at the saturation states across dimensions have been determined numerically [33,55]. These results reveal that saturated RSA packings are nonhyperuniform, even if the values of S(0) are relatively small (see Table I).

Poisson Cluster Process
The Poisson cluster process is an example of a strongly clustering point process with large density fluctuations on large length scales, i.e., with a large but finite value of S(0), and hence is a nonhyperuniform system that is far from being hyperuniform. The construction of the cluster process starts from a homogeneous Poisson point process of intensity ρ p [24]. Each point of the Poisson point process is the center of a cluster of points. The number of points in each cluster is independent and follows a Poisson distribution with mean value c. In our specific model, the positions of the points relative to the center of the cluster follows an isotropic Gaussian distribution with standard deviation r 0 , which can be regarded to be the characteristic length scale of a single cluster. This model is also known as a (modified) Thomas point process, which is an example of a Neyman-Scott process [36,57]. In the infinite-volume limit, the pair correlation function in R d is exactly given by [57]: Thus, the corresponding structure factor for any d is given by and hence such processes are nonhyperuniform and super-Poissonian with S(0) = 1 + c. To simulate the process, which is straightforward, we use periodic boundary conditions and the following parameters across the first three space dimensions: r 0 = 1 and unit number density ρ = ρ p c = 1 such that ρ p = 0.1 and c = 10. For such parameters, S(0) = 11 across dimensions (see Table I).

Hyperplanes Intersection Process
The hyperplanes intersection process (HIP) is hyperfluctuating [16], i.e., its number variances scales faster than the volume of the observation window and lim k→0 S(k) = ∞ [26,58]. This antihyperuniform and super-Poissonian point process is defined as the vertices (i.e., intersections) of a Poisson hyperplane process, that is, of randomly and independently distributed hyperplanes [36,59]. In the infinite-volume limit, the pair correlation function in R d for any d ≥ 2 is exactly given by [26]: where s is the specific surface of the hyperplane and ω d denotes the volume of a d-dimensional sphere of unit radius. The number density ρ is determined by the specific surface area s of the Poisson hyperplane process (which is the only parameter of the isotropic HIP): According to (77), because α = 1 for any d, the number variance has the large-R scaling σ 2 (R) ∼ R 2d−1 . Clearly, this process does not exist for d = 1. To simulate this process, we cannot employ periodic boundary conditions; rather, we circumscribe the cubic simulation box by a hypersphere and then generate intersecting hyperplanes that are Poisson distributed [58]. The orientation of the hyperplanes is uniformly distributed on the unit sphere and the distance of the hyperplanes to the center of the simulation box is uniformly distributed between zero and the radius of the circumsphere. The point process at unit number density is then simulated by computing all intersections of hyperplanes (within the circumsphere).

Hypercubic Lattice
Interestingly, the problem of determining number fluctuations in lattices has deep connections to number theory, including Gauss's circle problem [60] and its generalizations [16] as well as the Epstein zeta function [61], which is directly related to the minimization of the number variance [13,16,49]. All periodic point patterns in R d , including Bravais lattices, are hyperuniform of class I [13,16,49]. The hyperuniformity concept enables one to rank order lattices and other periodic point patterns according to the degree to which they suppress large-scale density fluctuations as defined by the number variance [13,16,49].
For the purposes of this investigation, it is sufficient to consider the higher-order fluctuations of the hypercubic lattice Z d is defined by where Z is the set of integers (. . . − 3, −2, −1, 0, 1, 2, 3 . . .) and x 1 , . . . , x d denote the components of a lattice vector.

Uniformly Randomized Lattice
It is well known that if the sites of a lattice are stochastically displaced by certain finite distances, the scattering intensity (structure factor) inherits the Bragg peaks (long-range order) of the original lattice, in addition to a diffuse contribution. It has recently been demonstrated that these Bragg peaks can be hidden in the scattering pattern for certain independent and identically distributed perturbations. We have referred to this protocol as the uniformly randomized lattice (URL) model [62]. The underlying long-range order can be "cloaked", under certain conditions, in the sense that it cannot be reconstructed from the pair-correlation function alone. Here we generate the URL model using the hypercubic lattice Z d and displace each lattice point by a random vector that is uniformly distributed in a rescaled fundamental cell of the lattice with aF ≡ [−a/2, a/2) d . The constant a controls the strength of perturbations. Counterintuitively, the long-range order suddenly disappears at certain discrete values of a and reemerges for stronger perturbations, as we will show. Here we cloak the Bragg peaks of Z d using the special value a = 1. Such cloaked URLs are hyperuniform such that S(k) ∼ k 2 in the limit k → 0 and hence are of class I [see Eq. (76)].

Stealthy Hyperuniform Process
Stealthy hyperuniform processes are defined by a structure factor that vanishes in a spherical region around the origin, i.e., S(k) = 0 for 0 < |k| ≤ K. Such point processes are hyperuniform of class I; see Eq. (76). A powerful procedure that enables one to generate high-fidelity stealthy hyperuniform point patterns is the collectivecoordinate optimization technique [63][64][65][66][67][68]. This optimization methodology involves finding the highly degenerate ground states of a class of bounded pair potentials with compact support in Fourier space, which are stealthy and hyperuniform by construction. The control parameter χ is a dimensionless measure of the ratio of constrained degrees of freedom (i.e., wave vectors contained within the cut-off wavenumber K) to the total degrees of freedom (approximately dN ) in such an optimization procedure. A point configuration with a small value of χ (relatively unconstrained) is disordered, and as χ increases, the short-range order increases within a disordered regime (χ < 1/2 for d = 2 and d = 3) [67]. For d = 1, stealthy hyperuniform states can be disordered for χ < 1/3 [68]. Here we use the collective-coordinate procedure to generate "entropically-favored" disordered stealthy point processes by first performing molecular dynamics simulations at sufficiently low temperatures and then minimizing the energy to obtain ground states with exquisite accuracy [68]. Importantly, stealthy states possess the bounded-hole property [46,48] and hence, as discussed in Sec. III C, P [N (R) = 0] = 0 for R > R max (χ), where R max (χ) is the radius of the largest hole in space dimension d, which depends on the control parameter χ.

V. GAUSSIAN DISTANCE METRIC
As noted in the Introduction, we are interested in ascertaining how large R must be such that P [N (R)] is well approximated by the normal distribution. There are several candidate "distance" metrics that we have considered that could be used to quantify proximity to the Gaussian distribution (beyond the skewness and excess kurtosis). One possible distance metric that we considered is the Kolmogorov-Smirnov test statistic [69]. While it is statistically robust, we found it be too insensitive for our purposes. The Kullback-Leibler divergence (also called the relative entropy) is a well-known measure of the difference between probability distributions [70]. However, it is not well-defined for comparing a discrete to a continuous distribution; see the SM for details.
A recent study that considered distance metrics between a pair of general functions that depend on ddimensional vectors was based on the integrated squared difference in R d , i.e., an L 2 distance metric [71]. These authors also found that the Kullback-Leibler distance was not useful for their purposes. These findings motivated us to consider distance metrics based on the squared difference between the Gaussian and number distributions.
We identify here two different contributions to the "distance" between a number distribution P [N (R)] and a Gaussian distribution: (1) deviations in the functional form of P [N (R)] and (2) the discreteness of N (R) that can only approximate the continuous Gaussian distribution. We found that the second contribution is essentially determined by the value of the number variance σ 2 (R), i.e., the contribution is smaller for larger values of σ 2 (R) for the following reason. Consider the standardized random variable [N (R)− N (R) ]/σ(R) whose discrete probability mass function is to approximate the continuous normal density. Then, the bin width of the probability mass function is given by 1/σ(R) and converges to zero for σ(R) → ∞.
Moreover, we found that the weight of the two contributions (relative to each other) strongly depends on the representation of the number distribution, e.g., via the characteristic function (Fourier representation) [69] or via direct space representations either in discrete or continuous forms. In fact, the choice of representation can virtually reverse the order of the distance metrics for our point processes. For example, a strong contribution (2) may result in a lower distance metric for the highly skew distribution of the HIP than for the Poisson process. Because contribution (2) of the discreteness of P [N (R)] is already essentially given by σ 2 (R), we here choose a representation that focuses on deviations in the functional form of P [N (R)] from that of a Gaussian random variable.
Therefore, we define an integer-valued random variable G(R), whose probability mass function P G(R) is proportional to a Gaussian distribution with the same first and second moment as our number distribution at radius R. We introduce a type of L 2 distance metric, denoted by l 2 (R), which is a Gaussian distance metric that employs the cumulative distribution function: where F G (n) is the cumulative distribution function of G(R), i.e., and F N (n) is the cumulative probability distribution of N (R), i.e., Note that the series of the squared differences in (81) is scaled by 1/σ(R) because for a Gaussian distribution the range of values for which F G ∈ [ε, 1 − ε] is proportional to σ(R). If a particular number distribution converges (sufficiently fast) to a Gaussian distribution, l 2 (R) will tend to zero. In summary, the Gaussian distance metric (81) is designed to be sensitive to small deviations in the distribution functions and at the same time robust against statistical fluctuations.

VI. SAMPLING OF MOMENTS AND NUMBER DISTRIBUTION
We sample number fluctuations within a spherical window of radius R for all models using a two-step procedure. First, we randomly place the observation window in the sample (using a uniform distribution for its center). Second, we determine the number of points N (R) within the observation window using periodic boundary conditions, except for the hyperplanes intersection process (HIP), as described in Sec. IV A 5. To reduce computational resources, we use the same centers for all radii that we consider. Relevant simulation parameters and properties for each of the models described above across dimensions are listed in Table I.
Thus, we empirically determine the probability mass function P [N (R)] and compute the mean value, variance, skewness, excess kurtosis. We determine the first four central moments using the unbiased estimators from Ref. [77]. Finally, we compute the l 2 distance metric, as described above. A source of systematic errors, which must be avoided for both the moments and distance measures, can arise when the number of observation windows N window per sample is too large. This effect can cause the distance metric l 2 (R) to artificially increase again for large radii. Therefore, we have used between 1 and 10 4 observation windows per sample depending on the system size and computational cost, so that the the systematic error is either avoided completely or smaller than the effects caused by statistical fluctuations. As a rule of thumb, we chose N window such that the volume fraction of the union of all observation windows of the largest window radius R max in one sample should not exceed 50%, i.e., 1 − exp [−N window v 1 (R max )/V ] < 0.5, where V is the volume of a single sample. Our estimators for l 2 (R) are statistically robust for values that are larger than the TABLE I. Simulation parameters all of the model point processes across the first three space dimensions considered in the work. Here, N is the average of point number inside a fundamental cell, Nc is the number of point patterns considered, N window is the number of observation windows per point pattern, and Rmax is the largest radius of an observation window. We have also indicated the values of the structure factors at the origin, S(0). In the cases of the equilibrium packings, they are obtained from the exact result for hard rods (d = 1) [72] and highly accurate approximations for the conditional nearest-neighbor function GP (r) in the limit r → ∞ [73,74]. Note that our computer simulation results are in excellent agreement with these analytical estimates of S(0) for equilibrium packings.  inverse of the total number n of observation windows, where n = N window × N c with N c being the number of configurations. For a finite number of samples, we empirically find that l 2 (R) typically cannot be smaller than approximately O(1/n). Therefore, we apply a datacut and only consider radii up to We apply the same datacut to the skewness and excess kurtosis [78]. We show all data without the datacut in the SM. All simulated data are available at a Zenodo repository [79].
Another obvious source of systematic errors can arise when the size of the sampling window is not much smaller than the size of the simulation box with a fixed number of particles, i.e., when canonical ensembles are employed. It is well-known that such finite-size effects lead to an underestimation of the local number variance σ 2 finite (R) when compared to its value in the thermodynamic limit. An empirical formula to estimate the first-order correction to the thermodynamic limit [17] shows that the error term is proportional S(0) (because there are larger fluctuations in the number of points per simulation box). Importantly, this implies that hyperuniform models, defined by lim k→0 S(k) = 0, are more robust against such finite-size effects.

VII. RESULTS
We describe results that we have obtained for the second, third and fourth cumulants, σ 2 (R), γ 1 (R) and γ 2 (R), as a function of the window radius R for all models across the first three space dimensions at unit number density (ρ = 1), as well as the corresponding full probability distributions P [N (R)] and the Gaussian distance metric l 2 (R). A testament to the high-precision of the data is the excellent agreement with rigorous bounds for these quantities for general cases as well as with exact results for the cases of packings for certain R reported in Sec. III (see the SM for details).
A. Cumulants σ 2 (R), γ1(R) and γ2(R) for the Models Across Dimensions Figure 2 shows the second, third and fourth cumulants as a function of R versus the window radius R for all models for d = 1, d = 2 and d = 3. The large-R asymptotic behavior of the variance is determined by the structure factor at the origin, S(0) [13]. As expected, the scaled number variance σ 2 (R)/v 1 (R) grows with R for the HIP process for two-dimensional (2D) and threedimensional (3D) cases. For all other nonhyperuniform models, σ 2 (R)/v 1 (R) asymptotes to a constant for large R in all dimensions. Of course, this scaled variance decreases with R for the three hyperuniform models (hy-percubic lattice, URL and stealthy systems).
For d ≥ 2, we have found via analyses given in the SM and immediately below that the skewness γ 1 (R) and excess kurtosis γ 2 (R) for our disordered hyperuniform systems vanish faster with increasing R than those for nonhyperuniform systems. Among all models studied, the quantities γ 1 (R) and γ 2 (R) vanish slowest for the antihyperuniform HIP models for d ≥ 2. The specific decay rates for all models are described below.
It is noteworthy that one-dimensional systems can present fluctuation anomalies not present in higher dimensions. For example, in the case of the integer lattice, the random variable N (R) can only take at most two values for any R, which of course is abnormally non-Gaussian. While the hypercubic lattice for d ≥ 2 never achieves a CLT (as discussed below), the variance is considerably broader than that for d = 1. Another anomalous category is class I hyperuniform systems, which have bounded variance for d = 1 [see Eq. (76)] and thus any such hyperuniform point process cannot obey a CLT because the standardized distribution cannot converge to a continuous distribution. This is clearly borne out by the distance metric plots, shown in Fig. 5 and Fig. S8 of the SM, for both the 1D disordered stealthy point process and integer lattice.
For the Poisson and super-Poissonian models (cluster and HIP) across dimensions, both the skewness and excess kurtosis decay monotonically; see Fig. 2. By contrast, the skewness and excess kurtosis oscillate about zero for the lattice, stealthy, equilibrium and RSA systems across all dimensions because all of them exhibit at least short-range order. It should not go unnoticed how the oscillations in γ 1 (R) and γ 2 (R) are related to short-or long-range order at the level of three-and fourbody correlation functions, g 3 and g 4 , as can be seen from the explicit formulas (22) and (23) for the skewness and excess kurtosis. In the instances of RSA and equilibrium packings, oscillations in γ 1 (R) and γ 2 (R) arise from strong short-range order exhibited by g 3 and g 4 .
The disordered hyperuniform systems that we consider, stealthy and URL point processes, have extraordinary number fluctuation behaviors. For 2D and 3D stealthy systems, both γ 1 (R) and γ 2 (R) strongly oscillate about zero; see Fig. 2. This suggests that at the level of g 3 , stealthy systems, counterintuitively, exhibit significant ordering on much larger length scales than the short-range order seen in the pair correlation function [67]. For 1D stealthy systems, γ 1 (R) and γ 2 (R) show even stronger oscillations than their higher-dimensional counterparts, indeed reflecting possible long-range order present in g 3 and g 4 . It is remarkable that the skewness and excess kurtosis can detect such anomalous long-range order that would not be expected based solely on the behavior of the pair correlation function. Another reason that supports the capacity of γ 1 (R) and γ 2 (R) to detect unusual long-range order is the cloaked URL model, which at the level of the pair correlation function would be considered to be highly disordered. Whereas γ 1 (R) for this model is a monotonically decreasing function of R, γ 2 (R) exhibits oscillations. This is entirely consistent with the fact that the periodicity of the underlying lattice is completely hidden at the level of the three-point correlation function but manifests itself for the first time in g 4 [62]. Figures 3 and 4 show representative configurations of all of the 2D models and their corresponding standardized number distributions. (We provide corresponding figures for all 1D and 3D models in the SM). Except for hypercubic lattices for any d and 1D hyperuniform systems of class I, all of the considered models across dimensions obey a CLT. It is noteworthy that for all models, except the hypercubic lattices, the number distribution functions P [N (R)] are unimodal, i.e., one with a single peak. (For small radii R 0.5, the number distributions are monotonically decreasing, which are still considered to be unimodal.) Recall that for all models, except the hypercubic lattice and the class I hyperuniform models in one dimension, the skewness and excess kurtosis tend to zero for large R for all dimensions. For well-behaved unimodal distributions, such a vanishing of both γ 1 (R) and γ 2 (R) indicates a tendency to a CLT. Indeed, this is confirmed by visual inspection of our corresponding evaluations of the number distributions for all dimensions; see Figs. 3 and 4. Our conclusions about the tendencies to CLTs for sufficiently large R are further confirmed by our evaluations of the Gaussian distance metric l 2 (R) for all models across dimensions, which are plotted in Fig. 5. For d = 2 and d = 3, one clearly sees that disordered hyperuniform point processes are better approximated by the normal distribution than their nonhyperuniform counterparts at a given large value of R. This is consistent with the visual inspections of the plots presented Figs. 3 and 4. From Fig. 5 and Fig. S8 of the SM, one can ascertain, for each model, a radius R 0 above which the distribution can be deemed to be approximately Gaussian, i.e., l 2 (R) is below some threshold for R > R 0 , implying that it is approximately determined by only the mean and the variance. Typically, we find that R 0 is orders of magnitude smaller for our hyperuniform and standard nonhyperuniform models than for our super-Poissonian models.

B. Number Distributions for the Models Across Dimensions
Our results clearly show that the hypercubic lattices across dimensions do not obey a CLT; see Figs. 2, 4 and 5. This is due to the fact that lattice points in general are "rigid" in the sense that fluctuations in N (R) are always stringently bounded from below and above for any value of R (as discussed in Sec. III C) due to their inherent long-range order. Hence, P [N (R)] is highly sensitive to the value of R, so that both γ 1 (R) and γ 2 (R) rapidly oscillate around zero, but the amplitudes of the oscillations do not vanish as R increases. It is already well-established that the local variance σ 2 (R) of lattices exhibits such rapid oscillations; see Ref. [13] and references therein. Because of the rapid oscillations in γ 1 (R), γ 2 (R), and l 2 (R) for hypercubic lattices, we represent the data by points instead of curves, which would require interpolation between the data points. Note that the corresponding distance metrics l 2 (R) are bounded from below for any value of R. These observations are consistent with the visual inspections of the number distribution shown in Fig. 4 for d = 2 and those for d = 1 and d = 3 given in the SM.

C. Gamma-Distribution Approximation for All Models Obeying a CLT
We have studied a variety of well-known closed-form probability distributions to determine those which best approximate the actual distributions for finite R for all of our models. Remarkably, we have ascertained that the gamma distribution provides a good approximation to the number distribution P [N (R)] for all models that obey a CLT across all dimensions for intermediate to large values of R. (Recall that the hypercubic lattices and 1D hyperuniform systems do not obey a CLT.) The gamma distribution is defined by where k and θ are shape and scale parameters, respectively, which are related to the mean and variance of P [N (R)] as follows: It immediately follows that the associated skewness and excess kurtosis can be expressed simply in terms of the mean and number variance, yielding   Figure 6 compares the gammma-distribution approximations to simulation data for a representative disordered hyperuniform model (stealthy), a standard nonhyperuniform model (RSA) and an antihyperuniform model (HIP) in two dimensions for R = 5. Visual inspection reveals the gamma distribution provides a good approximation in each case. Similar good agreement between the gamma-distribution approximation and simulation data was found for other models (not shown in Fig. 6) obeying a CLT across dimensions, as discussed in the SM.
Importantly, the approximation of P [N (R)] by a gamma distribution enables us to estimate the large-R scalings of γ 1 (R) and γ 2 (R) for all models across dimensions that obey a CLT. Specifically, we find γ 1 (R) ∼ R −1/2 and γ 2 (R) ∼ R −1 for the antihyperuniform HIP, γ 1 (R) ∼ R −d/2 and γ 2 (R) ∼ R −d for standard nonhype-runiform models, and γ 1 (R) ∼ R −(d+1)/2 and γ 2 (R) ∼ R −(d+1) for the hyperuniform models. Table II summarizes these scaling behaviors. The fact that the excess kurtosis decays to zero faster than the skewness for any model that obeys a CLT, whether hyperuniform or not, implies that the dominant asymptotic correction of the gamma distribution to a CLT for large R is determined by the skewness; see Appendix A for a proof. It is noteworthy that these predictions based on the gamma distribution are consistent with numerical findings for all nonhyperuniform systems and the antihyperuniform HIP using an independent method that employs certain running averages of γ 1 (R) and γ 2 (R); see the SM for details. Moreover, for 2D and 3D disordered hyperuniform (stealthy and URL) models, the predictions from the gamma-distribution approximations are also consistent  with the observed scalings of the skewness. Due to the strong oscillations in the excess kurtosis for 2D and 3D stealthy and URL processes described above, it was numerically difficult to definitively determine their scalings from the running-average method. It should not go unnoticed that the exact formulas for the skewness and excess kurtosis, Eqs. (34) and (35), respectively, of the Poisson distribution are consistent with the scalings predicted by the gamma-distribution approximation, lending additional validity to the latter. Furthermore, we have found that for all models across dimensions that obey a CLT, the running-average procedure yields scalings for l 2 (R) that agree with the corresponding gamma-distribution approximations (see SM), which are identical to the scalings for the skewness γ 1 (R). The proof that l 2 (R) has exactly the same scaling as γ 1 (R) for the gamma distribution is given in Appendix A. Furthermore, our approximations of P [N (R)] by gamma distributions are also consistent with CLTs, since they converge to the Gaussian distribution in the limit of k → ∞ (i.e, R → ∞), as described in Appendix A. According to Table II, convergence to a CLT is slowest for the HIP (proportional to R −1/2 for d ≥ 2), followed by the Poisson cluster process, and the Poisson process. RSA and equilibrium packings have the same scaling behaviors as the Poisson cluster and Poisson processes, but with smaller coefficients of proportionality. The convergence to a CLT is fastest for the disordered hyperuniform processes for d = 2 and d = 3 [80] D. Effect of Dimensionality on the Approach to a CLT For any particular d-dimensional model that eventually obeys a CLT, does the number distribution function P [N (R)] tend to Gaussian-like behavior faster as the space dimension d increases? As noted in the Introduction, this question can be answered by appealing to the decorrelation principle [29], which states uncon- The gamma distribution provides good approximations for all models that obey a CLT across the first three dimensions considered, including the ones not shown here.  strained correlations in disordered packings that exist in low dimensions vanish as d tends to infinity, and all higher-order correlation functions g n for n ≥ 3 may be expressed in terms of the number density ρ and pair correlation function g 2 . The decorrelation principle begins to manifest itself in low dimensions for disordered packings [31][32][33] but other disordered systems with strongly repulsively interacting particles, including fermionic [34] and Gaussian-core point processes [35]. We know that the number distribution function P [N (R)] generally involves certain integrals over all of the n-body correlation functions. Therefore, the decorrelation principle implies that for any d-dimensional model that decorrelates with d, P [N (R)] increasingly becomes Gaussian-like as d increases, since the first and second moments, determined by ρ and g 2 , dominate the distribution. By the same token, for any model that correlates with increasing d, P [N (R)] increasingly deviates from the normal distribution as d increases. We have verified these broad conclusions for the models studied in this article. In Fig. 7, we plot the Gaussian distance metric l 2 (R) for a representative disordered hyperuniform model (URL), a sub-Poisson nonhyperuniform model (RSA) and an antihyperuniform model (HIP) across the first three space dimensions, respectively. (Recall that the URL obeys a CLT for d ≥ 2.) As expected, we see that for both the URL and RSA models for a fixed value of R (for R > 1), l 2 (R) decreases with increasing d (increasingly tends to Gaussian-like behavior) because of decorrelation, while for the HIP, l 2 (R) increases with increasing d (moves away from Gaussian-like behavior) because it increasingly correlates with d. Importantly, when comparing fluctuations across dimensions, one must choose a meaningful length scale to make the window radius R dimensionless. A simple and good choice is ρ 1/d , which is proportional to the mean nearest-neighbor distance in space dimension d, and explains why the horizontal axes in each subfigure of Fig. 7 is ρ 1/d R.

VIII. CONCLUSIONS AND DISCUSSION
Via theoretical methods and high-precise simulation studies, we accurately quantified the skewness γ 1 (R), excess kurtosis γ 2 (R) and the number distribution P [N (R)] for eight different models of statistically homogeneous point processes in two and three dimensions: five nonhyperuniform models, one of which is anti-hyperuniform (HIP), and three hyperuniform models. Analogous models were also examined in one dimension, except for HIP, which is not defined in this dimension. We validated our simulation results for γ 1 (R), γ 2 (R) and P [N (R)] for all models by showing that they are in excellent agreement with rigorous bounds and exact results that we derived for the applicable ranges of R. For all disordered hyperuniform models, our explicit general formulas for γ 1 (R) and γ 2 (R) in terms of n-body information [Eqs. (22) and (23)] enable us to infer the existence of "hidden" type of long-range order that manifests itself for the first time at the three-body level or higher. Thus, the skewness and excess kurtosis have the capacity to detect anomalous long-range order that would not be expected based on the behavior of the pair correlation function alone.
We have introduced a novel Gaussian "distance" metric l 2 (R) to ascertain the proximity of P [N (R)] to the normal distribution for each model as a function of R. We have verified that l 2 (R) is a sensitive metric via numerical and theoretical methods. Since the distributions for all models (except the lattices) across dimensions are unimodal, the tendency to a CLT corresponds to the skewness and excess kurtosis simultaneously tending to zero. Almost all of the considered models across dimensions obey a CLT. We found that disordered hyperuniform point processes are better approximated by the normal distribution than their nonhyperuniform counterparts at a given large value of R. We proved that any 1D hyperuniform system of class I as well as the hypercubic lattice for any d cannot obey a CLT. Similarly, a general lattice in any dimension cannot obey a CLT.
It is noteworthy that we discovered that the gamma distribution provides a good approximation to the number distribution P [N (R)] for all models that obey a CLT across all dimensions for intermediate to large values of R, enabling us to estimate the large-R scalings of γ 1 (R), γ 2 (R) and l 2 (R). These predictions were corroborated by corresponding simulation results in almost all instances, as detailed in the SM. It is only in the special cases of the excess kurtosis for 2D and 3D stealthy and URL processes where the running-average method was not reliable enough to definitively determine their scalings due to strong oscillations, and thus this represents a simulation challenge for future research. Among all models, the convergence to a CLT is generally fastest for the disordered hyperuniform processes in two and three dimensions such that γ 1 (R) ∼ l 2 (R) ∼ R −(d+1)/2 and γ 2 (R) ∼ R −(d+1) for large R. The convergence to a CLT is slower for standard nonhyperuniform models such that γ 1 (R) ∼ l 2 (R) ∼ R −d/2 and γ 2 (R) ∼ R −d . Not surprisingly, the convergence to a CLT is slowest for the antihyperuniform HIP model such that γ 1 (R) ∼ l 2 (R) ∼ R −1/2 and γ 2 (R) ∼ R −1 . Using the decorrelation principle, we elucidated why any d-dimensional model that "decorrelates" or "correlates" with d corresponds to a P [N (R)] that increasingly moves toward or away from a CLT, respectively.
It has recently been reported that in a 2D absorbing phase of a lattice gas model, the density distribution at a hyperuniform length scale increasingly deviates from a Gaussian distribution on approach to the critical point [18], corresponding to a class III hyperuniform system [16]. These distributions are distinguished from those of our models in that they have heavy left tails (γ 1 (R) < 0). This implies that they cannot be approximated by gamma distributions, as in the majority of the models studied in the present paper across dimensions.
A related but distinctly different study from the one considered in the present work concerns fluctuations when the window is centered on a point of the point process [9,43]. It will be interesting in future work to carry out the analogous investigation of the corresponding skewness, excess kurtosis and number distribution for such particle-based quantities for hyperuniform and nonhyperuniform models.