Statistical properties of determinantal point processes in high-dimensional Euclidean spaces

Statistical properties of determinantal point processes in high-dimensional Euclidean spaces

Antonello Scardicchio ascardic@princeton.edu Department of Physics, Joseph Henry Laboratories, Princeton University, Princeton, New Jersey 08544, USA; and Princeton Center for Theoretical Science, Princeton University, Princeton, New Jersey 08544, USA    Chase E. Zachary czachary@princeton.edu Department of Chemistry, Princeton University, Princeton, New Jersey 08544, USA    Salvatore Torquato torquato@princeton.edu Department of Chemistry, Princeton University, Princeton, New Jersey 08544, USA;
Program in Applied and Computational Mathematics, Princeton University, Princeton, New Jersey 08544, USA;
Princeton Institute for the Science and Technology of Materials, Princeton University, Princeton, New Jersey 08544, USA;
Princeton Center for Theoretical Science, Princeton University, Princeton, New Jersey 08544, USA; and School of Natural Sciences, Institute for Advanced Study, Princeton, New Jersey 08540, USA
Abstract

The goal of this paper is to quantitatively describe some statistical properties of higher-dimensional determinantal point processes with a primary focus on the nearest-neighbor distribution functions. Toward this end, we express these functions as determinants of matrices and then extrapolate to . This formulation allows for a quick and accurate numerical evaluation of these quantities for point processes in Euclidean spaces of dimension . We also implement an algorithm due to Hough et. al. Hough et al. (2006) for generating configurations of determinantal point processes in arbitrary Euclidean spaces, and we utilize this algorithm in conjunction with the aforementioned numerical results to characterize the statistical properties of what we call the Fermi-sphere point process for to . This homogeneous, isotropic determinantal point process, discussed also in a companion paper Torquato et al. (2008), is the high-dimensional generalization of the distribution of eigenvalues on the unit circle of a random matrix from the circular unitary ensemble (CUE). In addition to the nearest-neighbor probability distribution, we are able to calculate Voronoi cells and nearest-neighbor extrema statistics for the Fermi-sphere point process and discuss these as the dimension is varied. The results in this paper accompany and complement analytical properties of higher-dimensional determinantal point processes developed in Torquato et al. (2008).

I Introduction

Stochastic point processes (PPs) arise in several different areas of physics and mathematics. For example, the classical statistical mechanics of an ensemble of interacting point particles is essentially the study of a random point process with the Gibbs measure providing the joint probability measure for an -tuple of vectors to be chosen. Moreover, some many-body problems in quantum mechanics, as we will see, can be regarded as stochastic point processes, where quantum fluctuations are the source of randomness. With regard to mathematical applications, it has been well-documented Mehta (2004) that the distribution of zeros of the Riemann zeta function on the critical line is well-represented by the distribution of eigenvalues of a random Hermitian matrix from the Gaussian unitary ensemble (GUE) or circular unitary ensemble (CUE) in the limit . Nevertheless, it remains an open problem to devise efficient Monte Carlo routines aimed at sampling these processes in a computationally efficient way.

In studies of the statistical mechanics of point-like particles one is usually interested in a handful of quantities such as -particle correlation functions, the distributions of the spacings of particles, or the distributions of the sizes of cavities. Although these statistics involve only a small number of particles, it is not simple to extract them from knowledge of the joint probability density . In general numerical techniques are required because analytical results are rare. It is then of paramount importance to study point processes for which analytic results exist for at least some fundamental quantities. The quintessential example of such a process is the so-called Poisson PP, which is generated by placing points throughout the domain with a uniform probability distribution. Such a process is completely uncorrelated and homogeneous, meaning each of the -particle distribution functions is equal to , where is the number density for the process. Configurations of points generated from this process are equivalent to classical systems of noninteracting particles or fully-penetrable spheres Torquato (2002), and almost all statistical descriptors may be evaluated analytically.

One nontrivial example of a family of processes which has been extensively studied is the class of determinantal PPs, introduced in 1975 by Macchi Macchi (1975) with reference to fermionic statistics. Since their introduction, determinantal point processes have found applications in diverse contexts, including random matrix theory (RMT), number theory, and physics (for a recent review, see Soshnikov (2000)). However, most progress has been possible in the case of point processes on the line and in the plane, where direct connections can be made with RMT Mehta (2004) and completely integrable systems Jimbo et al. (1980).

Similar connections have not yet been found, to the best of our knowledge, for higher dimensional determinantal point processes, and numerical and analytical results in dimension are missing altogether. In this paper and its companion Torquato et al. (2008), we provide a generalization of these point processes to higher dimensions which we call Fermi-sphere point processes. While in Torquato et al. (2008) we have studied, mainly by way of exact analyses, statistical descriptors such as -particle probability densities and nearest-neighbor functions for these point processes here we base most of our analysis on an efficient algorithm Hough et al. (2006) for generating configurations from arbitrary determinantal point processes, and are therefore able to study other particle and void statistics related to nearest-neighbor distributions and Voronoi cells.

In particular, after presenting in detail our implementation of an algorithm Hough et al. (2006) to generate configurations of homogenous, isotropic determinantal point processes, we study several statistical quantities thereof, including Voronoi cells statistics and distributions of minimum and maximum nearest neighbor distances (for which no analytical results exist). Additionally, the large- behavior of the nearest-neighbor functions is computationally explored. We provide substantial evidence that the conditional probabilities and , defined below, are asymptotically linear, and we give estimates for their slopes as a function of dimension between one and four.

The plan of the paper is as follows. Section II provides a brief review of determinantal point processes and defines the statistical quantities used to characterize these systems. Of particular importance is the formulation of the probability distribution functions governing nearest-neighbor statistics as determinants of matrices; the results are easily evaluated numerically. The terminology we develop is then applied to the statistical properties of known one- and two-dimensional determinantal point processes in Section III. Section IV discusses the implementation of an algorithm for generating determinantal point processes in any dimension , and we combine the results from this algorithm and the numerics of Section II to characterize the so-called Fermi sphere point process for and . In Section V we provide an example of determinantal point-process on a curved space (a 2-sphere) and our conclusions are collected in Section VI.

Ii Formalism of determinantal point processes

ii.1 Definitions: -particle correlation functions

Consider point particles in a subset of -dimensional Euclidean space . It is convenient to introduce the Hilbert space structure given by square integrable functions on ; we will adopt Dirac’s bra-ket notation for these functions. Unless otherwise specified, all integrals are intended to extend over . A determinantal point process can be defined as a stochastic point process such that the joint probability distribution of points is given as a determinant of a positive, bounded operator of rank :

(1)

where is the kernel of . In this paper, we focus on the simple case in which the nonzero eigenvalues of are all 1; the more general case can be treated with minor changes Hough et al. (2006). We can write down the spectral decomposition of as:

(2)

where are the eigenvectors of the operator . The reason for the superscript on the basis vectors will be clarified momentarily. The correct normalization of the point process is obtained easily since Mehta (2004):

(3)

where the last determinant is to be interpreted as the product of the non-zero eigenvalues of the operator . Since these eigenvalues are all unity we obtain , which yields:

(4)

Notice that in terms of the basis we can also write:

(5)

An easy proof is obtained by considering the square matrix . Then,

(6)

which is the same as (1).

Determinantal point processes are peculiar in that one can actually write all the -particle distribution functions explicitly. The -particle probability density, denoted by , is proportional to the probability density of finding the first particles in volume elements around the given positions , irrespective of the remaining particles. For a general determinantal point process this function takes the form:

(7)

In particular, the single-particle probability density is:

(8)

This function is proportional to the probability density of finding a particle at , also known as the intensity of the point process. One can see that the normalization is:

(9)

For translationally-invariant processes , independent of . We remark in passing that for a finite system translational invariance is defined in the sense of averaging the location of the origin over with periodic boundary conditions enforced.

It is also possible to write the two-particle probability density explicitly:

(10)

which has the following normalization:

(11)

In general the normalization for is given by , or the number of ways of choosing an ordered subset of points from a population of size . For a translationally-invariant and completely uncorrelated point process (10) simplifies according to .

We also introduce the -particle correlation functions , which are defined by:

(12)

Since for a completely uncorrelated point process, it follows that deviations of from unity provide a measure of the correlations between points in a point process. Of particular interest is the pair correlation function, which for a translationally-invariant point process of intensity can be written as:

(13)

Closely related to the pair correlation function is the total correlation function, denoted by ; it is derived from via the equation:

(14)

where the second equality applies for all determinantal point processes by (13). Since as () for translationally invariant systems without long-range order, it follows that in this limit, meaning that is generally an function, and its Fourier transform is well-defined.

Determinantal point processes are self-similar; integration of the -particle probability distribution with respect to a point gives back the same functional form 111One could think in terms of effective interactions and renormalization group. The determinantal form of the probabilities then is a fixed point of the renormalization operation of integrating out one or more particles.. This property is desirable since it considerably simplifies the computation of many quantities. However, we note that even complete knowledge of all the -particle probability distributions is not sufficient in practice to generate point processes from the given probability . This notoriously difficult issue is known as the reconstruction problem in statistical mechanics Torquato and Stillinger (2003, 2002a); Kuna et al. (2007); Uche et al. (2006). When in Section IV we discuss an explicit constructive algorithm to generate realizations of a given determinantal process, the reader should keep in mind that the ability to write down all the -particle correlation functions is not the reason why there exists such a constructive algorithm.

ii.2 Exact results for some statistical quantities

We have seen that the determinantal form of the probability density function allows us to write down all -particle correlation functions in a quick and simple manner. However, we can also express more interesting functions, such as the probability of having an empty region or the expected number of points in a given region, as properly constructed determinants of the operator . This property has been used in random matrix theory to find the exact gap distribution of eigenvalues on the line in terms of solutions of a nonlinear differential equation Tracy and Widom (1994a). The relevant formula is a special case of the result Soshnikov (2000) that the generating function of the distribution of the number points in the region is:

(15)

where is the characteristic function of , is the identity operator, and . We will also denote . Therefore, the probability that the region is empty is obtained by taking the limit in the previous formula. The result is:

(16)

Equation (16) may be written more explicitly. Consider the eigenvalues of . By the definition of the determinant, equation (16) takes the form:

(17)

where the product is over the non-zero eigenvalues of only (of which there are , the number of particles). First notice that for the non-zero we have , where are the eigenvalues of . In fact one can show that the traces of all the powers of these two operators are the same using , and the cyclic property of the trace operation. This condition is sufficient for finite, and the limit can be taken afterward. The operator can now be written in a basis as the matrix:

(18)

and the determinant in (16) as:

(19)

We will be using this formula often in the following analysis. The probability has a unique role in the study of various point processes Torquato (2002), in particular when , a ball of radius (for translationally-invariant processes the position of the center of the ball is immaterial). In this context, is called the void exclusion probability Torquato et al. (1990); Torquato and Lee (1990); Lu and Torquato (1992); Torquato (2002), and we will adopt this name and notation in this paper (in Torquato et al. (2008) we have studied this quantity in an appropriate scaling limit, when ).

However, there are statistical quantities of great importance which cannot be found with the above formalism. For example, one can examine the distribution of the maximum or minimum nearest neighbor distances in a determinantal point process, or the “extremum statistics,” and these quantities cannot be found easily by the above means. One could also explore the distribution of the Voronoi cell statistics or the percolation threshold for the PP. To determine these quantities we will have to rely on an explicit realization of a determinantal point process. The existence and the analysis of an algorithm to perform this task is a central topic of this paper.

We introduce now some quantities which characterize a PP Torquato et al. (1990); Torquato and Lee (1990); Lu and Torquato (1992); Torquato (2002). We start with the above expression for the probability of finding a spherical cavity of radius in the point process. Analogously, one can define the probability of finding a spherical cavity of radius centered on a point of the process, which we denote as . can be found in connection with using the following construction. Consider the probability of finding no points in the spherical shell of inner radius and outer radius , which we call . This function can be obtained by either of the previous formulas (16) or (19). It is clear that . It is also true that for sufficiently small the probability of having two or more points in the sphere of radius is negligibly small compared to the probability of having one particle. Hence, the probability of finding no particles in the spherical shell conditioned on the presence of one point in a sphere of radius and volume is:

(20)

and by taking the limit of this expression we find that:

(21)

That can be seen from the following argument: set . Then, because the region is infinitesimal and hence empty with probability 1, and since for sufficiently small we have at most one point in the region. One line of algebra provides the result.

Using this expression, we can derive an interesting and practical result for . First, notice that contains the matrix defined by (18), which when becomes:

(22)

Moreover, if we assume that is invertible, we can see that to first order in :

(23)

From (23) we find the final result:

(24)

where . Notice that for , we have , and as expected.

These two primary functions can be used to define four other quantities of interest. Two are density functions:

(25)
(26)

which can be interpreted as the probability densities of finding the closest particle at distance from a random point of the space or another random point of the process, respectively. The other two functions are conditional probabilities:

(27)
(28)

which give the density of points around a spherical cavity centered respectively on a random point of the space or on a random point of the process. We note that is the surface area of the -dimensional sphere of radius . We will study the behavior of these functions for some determinantal PPs in Sections III and IV of this paper.

From the definitions in (25)-(28) in conjuction with (19) and (24), it is possible to express and as numerically-solvable operations on matrices. The results are:

(29)
(30)
(31)
(32)

The form (which serves as a definition of ) in (32) is of particular interest. If the correction term for all , positivity and monotonicity of (which must be proven independently) are then sufficient to ensure that for appropriately large , in scaling. Although we have been unable to develop analytic results for the large- behavior of , numerical results, which are provided later (see Fig. 10), suggest that and monotonically as for , and constant for . As both behaviors are subdominant with respect to the linear growth of , we expect that and possess the same linear slope for sufficiently large .

An important point to address is the convergence of the results from (19) in the limit . We expect that the calculations for finite but large provide an increasingly sharp approximation to the results from the limit. Figure 1 presents the calculation of for a few values of with ; it is clear that the numerical calculations quickly approach a fixed function for , and it is this function which we accept as the correct large- limit. The results for higher-dimensional processes are similar, and we will assume that this convergence property holds throughout the remainder of the paper.

Figure 1: Convergence of the numerical results using (19) for with respect to increasing matrix size .

ii.3 Hyperuniformity of point processes

Of particular significance in understanding the properties of determinantal point processes is the notion of hyperuniformity, also known as superhomogeneity. A hyperuniform point pattern is a system of points such that the variance of the number of points in a spherical window of radius obeys:

(33)

for large Torquato and Stillinger (2003). This condition in turn implies that the structure factor has the following small- behavior:

(34)

meaning that hyperuniform point patterns do not possess infinite-wavelength number fluctuations Torquato and Stillinger (2003). Examples of hyperuniform systems include all periodic point processes Torquato and Stillinger (2003), certain aperiodic point processes Torquato and Stillinger (2003); Gabrielli et al. (2003), one-component plasmas Torquato and Stillinger (2003); Gabrielli et al. (2003), point processes associated with a wide class of tilings of space Gabrielli and Torquato (2004); Gabrielli et al. (2008), and certain disordered sphere packings Torquato and Stillinger (2006a); Torquato et al. (2008); Torquato and Stillinger (2002a, b). It has also been shown Torquato et al. (2008) that the Fermi sphere determinantal point process, described below, is hyperuniform.

The condition in (34) suggests that for general translationally invariant nonperiodic systems:

(35)

for some . However, hyperuniform determinantal point processes may only exhibit certain scaling exponents . One can see for a determinantal point process that:

(36)

where denotes the Fourier transform. Equation (36) therefore suggests that:

(37)

Taking the inverse Fourier transform of (37) gives the following large- scaling of :

(38)

The negative coefficient and the negative argument of the gamma function in (38) are crucially important. Since for all , it must be true , and this condition restricts the possible values of the scaling exponent . Namely, the behavior of the gamma function requires that fall into one of the intervals , , , and so forth. We remark that the integer-valued endpoints of these intervals are indeed valid choices for and imply that for sufficiently large . These values of are therefore types of “limiting values” that overcome the otherwise dominant asymptotic scaling of . We provide an example of a determinantal point process with the critical scaling in Section III.C; the resulting large- behavior for is seen to be exponential.

Iii Properties of known determinantal point processes

iii.1 One-dimensional processes

By far, the most widely studied examples of determinantal point processes are in one-dimension. In fact, the connection to (RMT) led others to explore the statistical properties of these systems even prior to the formal introduction of determinantal point processes. To make this connection explicit, consider an random Gaussian Hermitian matrix, i.e., a matrix whose elements are independent random numbers distributed according to a normal distribution. This class of matrices defines the Gaussian unitary ensemble or GUE. It is possible to see Mehta (2004) that the distribution induced on the eigenvalues of these random matrices is:

(39)

where is an appropriate normalization constant. By a standard identity for the Vandermonde determinant:

(40)

and by combining the rows of the matrix appropriately we find:

(41)

where the functions are the Hermite orthogonal polynomials normalized such that the coefficient of the highest power of is unity. Taking into account the weight , we can write in agreement with (5):

(42)

where the orthonormal basis set is:

(43)

is a normalization factor. Therefore, this distribution is equivalent to the the one induced by a system of non-interacting, spinless fermions in a harmonic potential. We note without proof that the other canonical random matrix ensembles (GOE and GSE) can also be expressed as determinantal point processes by introducing an internal vector index for the basis functions Mehta (2004); Soshnikov (2000).

Another prominent example of a determinantal point process is given by the unitary matrices distributed according to the invariant Haar measure; the resulting class is termed the circular unitary ensemble or CUE Mezzadri (2006). The eigenvalues of these matrices can be written in the form with ; they are distributed according to (5) with the basis:

(44)

Notice that the eigenvalues represent the positions of free fermions on a circle, where the Fermi sphere has been filled continuously from momentum 0 to .

Another possible one-dimensional process is obtained by changing the exponent in (39) to an arbitrary polynomial. This generalization has interesting connections to the combinatorics of Feynman diagrams and to random polygonizations of surfaces Di Francesco et al. (1995). For other examples of one-dimensional determinantal point processes we refer the reader to Soshnikov (2000).

iii.2 Exact results in one dimension

For historical reasons, the most-studied descriptor of determinantal point processes is the gap distribution function, which represents the probability density of finding a chord of length separating two points in the system for ; we denote this function by . For canonical ensembles of random matrices exact solutions for have been written in terms of solutions of well-known nonlinear differential equations Mehta (2004). We start with the following observation: after an appropriate rescaling of the eigenvalues, the gap distribution of eigenvalues of a random matrix is a universal function, depending only on the “nature” of the ensemble (unitary, orthogonal or symplectic) which defines the small- behavior of . For example, the two ensembles GUE and CUE defined above will have the same gap distribution in the limit . In the case of the GUE the limit is taken for the eigenvalues:

(45)

where is in the “bulk” of the distribution ( for large). One can prove that all the eigenvalues of a large random matrix will fall in an interval of size with probability 1 in the large- limit. After this rescaling, the kernel converges to the “sine kernel” in the large- limit Tracy and Widom (1994a, b):

(46)

From this result one can find the -particle correlation functions. In particular one finds for :

(47)

Applying this procedure to the CUE leads to the very same kernel; for a wider class of examples relevant to physics see Brézin and Zee (1993). Convergence of the kernel implies weak convergence of all the -particle correlation functions to universal distributions. These distributions are defined by the sine kernel, one of a small family of kernels which appear to be universal Tracy and Widom (1994a, b) in controlling large- limits of various statistical quantities of apparently different distributions. The study of the analytic properties of the kernels in this family yields a complete solution for the Janossy probabilities and edge distributions in one-dimensional systems.

Once the limiting kernel is identified, a solution for the gap distribution still requires a detailed mathematical analysis Tracy and Widom (1994a). An approximate form for , known as Wigner’s surmise, was suggested by Wigner in 1951:

(48)

and it is an extremely good fit for numerical data. However, our primary focus in this work is on the asymptotic behavior of the conditional probability , and we therefore look for an exact solution for this function. First, we note without proof Forrester and Witte (2000) that for may be expressed in terms of a Painleve V transcendent. Namely, let be a solution of the nonlinear equation:

(49)

subject to the boundary condition:

(50)

as . We may then write in the form:

(51)

We recall that may also be expressed in terms of via the following relation:

(52)

By making a change of variables and comparing (51) and (52), we conclude that:

(53)

Equation (53) allows us to develop small- and large- expansions of in terms of the equivalent expansions for .

To describe the small- behavior of , we substitute an expansion of the form:

(54)

into (49) and solve order-by-order for the coefficients . Upon converting the solution to a result for using (53), we obtain:

(55)

The derivation of the large- expansion is similar. We choose an expansion of the form:

(56)

and substitute this equation into (49). After converting the result to an asymptotic series for with (53), we obtain:

(57)

By looking at Figure 2, one can see that the expansions are quite good for the ranges in where they are valid. Equations (49), (III.2), and (57) constitute the solution to our problem. Although it is natural to ask if there is a corresponding nonlinear differential equation that characterizes in higher dimensions, we are not aware of any work in this direction, and this issue remains an open problem.

Figure 2: Comparison of the exact form of for the determinantal point process with the small- and large- expansions in (III.2) and (57).

iii.3 Two-dimensional processes

There are a few examples of determinantal point processes in two dimensions. The seminal example is provided by the complex eigenvalues of random non-Hermitian matrices Ginibre (1965); Forrester (1998). The kernel of such a determinantal point process is given by:

(58)

where is the rank of the matrix and . Incidentally, (58) can be related to the distribution of polarized electrons in a perpendicular magnetic field, filling the lowest Landau levels. In the limit (58) becomes

(59)

which is a homogeneous and isotropic process () in . It is instructive to examine the pair correlation function, which after some algebra can be written as:

(60)

From this expression one finds that the correlation between two points decays like a Gaussian with respect the distance separating the points. Letting , we may write the associated structure factor of the system as:

(61)

which has the following small- behavior:

(62)

We see that the determinantal point process generated by the Ginibre ensemble is hyperuniform with an exponential scaling for small , corresponding to an endpoint of one of the “allowed” intervals for determinantal PPs; the large- behavior of the kernel is exponential ().

Other ensembles of two-dimensional determinantal point processes can be found in simple systems. For example, the zeros of an analytic Gaussian random function also form a determinantal point process on the open unit disk Peres and Virág (2000); Leboeuf (1999). The limiting kernel governing these zeros is called the Bargmann kernel:

(63)

and is inherently different from (59).

Iv An algorithm for generating determinantal point processes

We are able to write down an algorithm, which we call the HKPV algorithm, after Hough et al. (2006), to generate determinantal point processes due to the geometric interpretation of the determinant in as the volume of the simplex built with the vectors . In the original paper Hough et al. (2006) this algorithm is sketched and then proved to produce the correct distribution function . The algorithm is extremely powerful and versatile and we believe it is important to provide as many details as possible about it and its implementantion (which has not been done before, to our knowledge). Therefore we dedicate the present section to provide a complete description of the HKPV algorithm and enough details (with some tricks) for its efficient implementation.

Set , the kernel of the determinantal point process. Pick a point distributed with probability:

(64)

With this point build the new operator , defined by:

(65)

This operator has with probability 1 a single nonzero eigenvalue and null eigenvalues. When expressed as a matrix in the basis , takes the form:

(66)

Consider the null eigenvectors of ; we will denote them as and call the only eigenvector with a nonzero eigenvalue. The null eigenvectors can be found easily by means of a fast routine based on singular value decomposition (SVD), but we will see one that can proceed without it.

Next, build the new operator :

(67)

To simplify the computation, notice that by completeness of the basis in the eigenspace of :

(68)

and since is the only eigenvector orthogonal to the null space:

(69)

From this equation we conclude:

(70)

Once is obtained, we repeat the procedure with , generating the point from the probability distribution:

(71)

and the operators . As the number of iterations increases, we constantly reduce the rank of the operators by one: , , etc. Therefore, after we have placed the last point , we are left with the an operator of rank 0, and the algorithm stops. Reference Hough et al. (2006) shows that the -tuples are distributed according to the distribution (1).

The whole procedure requires steps for every realization, which is equal to the number of function evaluations necessary to create the matrices . Therefore, the algorithm is computationally quite light. The only subroutine that requires some work is the extraction of the random points from the probability distributions . For one can use a numerically-implemented inverse CDF technique Devroye (1986), and the computational cost of this procedure is independent of . For if the distributions are not very peaked, a rejection algorithm is sufficient. The rejection algorithm works by sampling points from a uniform distribution on the domain. A tolerance value near the maximum of the probability density of the point process is set, and the point is accepted if a uniform random number chosen between 0 and the tolerance value is less than the probability density at that point. Otherwise, the point is rejected, and the process repeats. Unfortunately, it is difficult to estimate the computational cost of this algorithm as a function of the number of particles .

iv.1 Numerical results in one dimension

We have implemented the algorithm described above to study a determinantal point process on the circle , where are the orthonormal functions with . This ensemble, as mentioned above, is equivalent to the one generated by the eigenvalues of unitary random matrices chosen according to the Haar measure. Eigenvalues of matrices from the CUE can be generated easily by means of a fast SVD algorithm Mezzadri (2006); however, plenty of exact results exist. Therefore, we study this one-dimensional determinantal process as a test both for the performance of our algorithm and for the convergence of the results to the limit.

We have implemented the algorithm in both Python and C++, noticing little difference in the speed of execution, and run it on a regular desktop computer. As mentioned above, the algorithm runs polynomially in with the sampling of the distribution limiting the computational speed. One point, however, which requires attention is the loss of precision of the computation. Due to the fixed precision of the computer calculations, the matrix ceases to have exactly integral trace, diminishing the reliability of the results. Typically, one observes deviations in the 5th decimal place after particles have been placed. We have devised an ‘error correction’ procedure in which the numerical matrix is projected onto the closest Hermitian matrix which has eigenvalues 1 or 0 only. This projection corrects for a great part, but not all, of the error; however, the algorithm is slowed by this modification. The number of particles in each configuration can therefore be pushed to , regardless of the dimension. We have been able to generate between 75000 and 100000 configurations of points in each dimension. In general, the error-correction procedure generates more reliable statistics for a given value of compared to the uncorrected algorithm, and we therefore expect that any residual error not captured by the matrix projection is minimal.

Figure 3: Comparison of the exact expression (75) for with the results from the HKPV Algorithm for The results from the simulation are obtained using 75000 configurations of 45 particles.

As a preliminary check for the HKPV Algorithm, we have calculated the pair correlation function and compared the results to the exact expression in (75) below. The comparison is quite favorable and suggests that the point configurations are being generated correctly by the algorithm. We mention a few characteristics of which arise from the determinantal nature of the point process. First, the system is strongly correlated for a significant range in , and as . This correlation hole McWeeney (1960); Slater (1951); Boyd and Coulson (1974); Schwabl (2005) is indicative of a strong effective repulsion in the system, especially for small point separations. In other words, the points tend to remain relatively separated from each other as they are distributed through space. Second, for all , meaning that it is always negatively correlated; again, this quality is indicative of repulsive point processes, which are characterized by a reduction of the probability density near each of the coordination shells in the system. We show in a separate paper Torquato et al. (2008) that at fixed number density, approaches an effective pair correlation function as , suggesting that the points achieve an increasingly strong effective hard core as the dimension of the system increases. At fixed mean nearest-neighbor separation this observation implies that for all as (as for any ), implying that the points become completely uncorrelated in this limit. We will show momentarily that the latter limit is difficult to interpret due to the dimensional dependence of the density .

Figure 4 presents the results for the gap distribution function using both the HKPV Algorithm and a numerical calculation based on the determinant in (19). As with the calculation of , the comparison between the numerical results and the simulation is favorable. This curve, as expected, has the same form as the one reported in the random matrix literature Mehta (2004) and scales with as . We stress, however, that this function represents the distribution of gaps between points on the line and does not discriminate between gaps to the left and to the right of a point. The random matrix literature oftentimes describes this quantity as a “nearest-neighbor” distribution, which it is not. As mentioned in the discussion following (25) and (26), the void and particle nearest-neighbor distribution functions are given by the functions and , respectively, and require that distance measurements be made both to the left and to the right of a point; the numerical and simulation results for these functions are also given in Figure 4.

Figure 4: Comparison of numerical and simulation results with for left: the gap distribution function , center: and right: .

The function is clearly different from . has a similar shape to the gap distribution function; however, peaks more sharply around while has a less intense peak near . This observation is justified from a numerical standpoint since point separation measurements are made in both directions from a given reference point with only the minimum separation contributing to the final histogram of . In contrast, every gap in the point process is used for constructing the histogram of . As a result, we expect the first moment of to be less than that of , and this result is exactly what we observe in Figure 4.

Figure 5: Numerical results using (19) for and with . Also included are representative simulation results and estimated errors from the HKPV Algorithm under the same conditions as in Figure 3.

The form of may at first seem confusing in the context of our discussion above concerning the inherent repulsion of the determinantal point process. Unlike and , the void nearest-neighbor function has a nonzero value at the origin and is monotonically decreasing with respect to . To understand this behavior, it is useful to examine the behavior of the corresponding and functions, which are plotted in Figure 5. We recall from (27) and (28) that and are related to conditional probabilities which describe, given a region of radius empty of points (other than at the center for ), the probability of finding the nearest-neighbor point in a spherical shell of volume , where is the surface area of a -dimensional sphere of radius . Of particular relevance to the behavior of is the fact that and for . Therefore, the dominant factor controlling the small- behavior of is the spherical surface area Torquato et al. (2008). Since is nonzero for , it follows from (27) that is nonzero in contrast to .

The behavior of both and is of particular interest in this paper. We conjecture that both functions are linear for sufficiently large in any dimension. We show elsewhere Torquato et al. (2008) that as , and , where is a dimensionally-dependent constant (for this is evident in Figure 5). Additionally, we believe that and obtain the same slope in the large- limit, and we will provide further commentary on this notion momentarily (see Fig.10). It is clear from Figure 5 that the results from the simulations are in agreement with the numerical results for a wide but limited range of and they begin to deviate respectably for sufficiently large. This is due to the fast decay of both and to zero, therefore giving very small statistics (and a large degree of uncertainty) at these values of . This said, the numerical results are clear and provide strong support for our claims above.

iv.2 Fermi sphere determinantal point process for

iv.2.1 Definition of the Fermi sphere point process

Here we study the determinantal point process of free fermions on a torus, filling a Fermi sphere. A detailed description of this process in any dimension may be found in an accompanying paper Torquato et al. (2008). We consider this example because it is the straightforward generalization of the one-dimensional CUE process described above. However, sampling of this ensemble cannot be accomplished with methods other than the algorithm introduced above; this limitation is in contrast to the two examples from Section III.C, where the ensemble may be generated from zeros of appropriate random complex functions. Nevertheless, it is difficult to construct another procedure that can be generalized to higher dimensions since zeros of complex functions and random matrices are naturally constrained to .

We consider the determinantal point process obtained by “filling the Fermi sphere” in a -dimensional torus, i.e., ; our choice of the box size is for convenience and without loss of generality. We therefore consider all functions of the form:

(72)

with

(73)

where is implicitly defined by the total number of states contained in the reciprocal-space sphere. This process is translationally-invariant for any , both finite or infinite, and isotropic in the limit ; it possesses the symmetry group of the boundary of the set (73), a dihedral group which approximates very well for sufficiently large. The pair correlation function can be easily calculated for any , and it is well-defined in the thermodynamic limit:

(74)

where and are -dimensional real vectors, and the sums extend over the set (73), which contains points. In the limit the sums become integrals over a sphere of radius , where is the number density. The resulting pair correlation function is given by:

(75)

where is the Bessel function of order (cf. Torquato et al. (2008)). This pair correlation function is clearly different from (60) for ; the two are therefore not equivalent, even in the thermodynamic limit. One can also find the limiting kernel:

(76)

which is also different from (59) and (63) for .222Different fillings of the spectrum give raise to different families of determinantal point processes. In Torquato et al. (2008) we present another example of determinantal point process where states with momenta where or are filled. We called these Fermi-shell point processes.

Figure 6 shows configurations of points generated for the and Fermi sphere point process alongside corresponding configurations for the Poisson point process in these dimensions. The repulsive nature of the determinantal point process is immediately apparent from these figures; note especially that the Fermi sphere point process discourages clustering of the points in space. In contrast, clustering is not prohibited in the Poisson point process, and small two- and three-particle clusters are easily identified. Of particular interest is that the Fermi sphere point process distributes the points more evenly through space due to the effective repulsion in the system. This characteristic reflects the hyperuniformity of the point process Torquato and Stillinger (2003), and we will have more to say about this property momentarily.

Figure 6: Upper: A configuration of points distributed according to left: a Fermi sphere determinantal point process and right: a Poisson point process. Lower: A configuration of