Realizable Hyperuniform and Nonhyperuniform Particle Configurations with Targeted Spectral Functions via Effective Pair Interactions

At first glance, the aforementioned Fourier-based collective-coordinate optimization procedure may seem to be ideally suited to construct possibly realizable hyperuniform configurations, since it enables one to obtain configurations with desired structure factors for wave vectors around the origin with very high accuracy Uche et al. (2006b); Batten et al. (2008); Zachary and Torquato (2011); Zhang et al. (2016). One begins with a single classical configuration of $N$ particles with positions $\mathbf{r}^{N}\equiv{\bf r}_{1},\ldots,{\bf r}_{N}$ in a fundamental cell under periodic boundary conditions. Here, the structure factor of a single configuration, $\mathcal{S}(\mathbf{k})$, is constrained to be equal to a target function $\mathcal{S}_{0}(\mathbf{k})$ for $\mathbf{k}$ in a certain finite set $\mathbb{K}$. These constraints are enforced by minimizing a fictitious potential energy $\Phi(\mathbf{r}^{N})$, defined to be the square of the difference between $\mathcal{S}({\bf k})$ and $\mathcal{S}_{0}({\bf k})$:

$$\Phi(\mathbf{r}^{N})=\sum_{\mathbf{k}\in\mathbb{K}}[\mathcal{S}(\mathbf{k})-\mathcal{S}_{0}(\mathbf{k})]^{2},$$

(15)

where, for a single configuration, the structure factor at a non-zero $\mathbf{k}$ vector is given by

$$\mathcal{S}(\mathbf{k})=\frac{1}{N}\left|{\tilde{\rho}}(\mathbf{k})\right|^{2}\mbox{ ($\mathbf{k}\neq\mathbf{0}$)},$$

(16)

and

$${\tilde{\rho}}(\mathbf{k})=\sum_{j=1}^{N}\exp(-i\mathbf{k}\cdot\mathbf{r}_{j})$$

(17)

is the complex collective density variable. Throughout the paper, we will use $\mathcal{S}$ to denote single-configuration structure factors and $S$ to denote ensemble-averaged structure factors. It was shown that the potential energy given by (15) is equivalent to a certain combination of a long-ranged two-, three-, and four-body interactions Uche et al. (2006b). Therefore, constraining $\mathcal{S}(\bf k)$ to a target function $\mathcal{S}_{0}(\bf k)$ using this method is equivalent to finding a single ground-state configuration with these interactions. One calculates $\mathcal{S}(\mathbf{k})$ from (16) rather than (2) not only because it applies only for a single finite-size configuration (not an ensemble), but (16) allows $\mathbb{K}$ to include $\mathbf{k}$ vectors very close to the origin.

However, this standard collective-coordinate method cannot be used for the realizability problem. It suffers from numerical difficulty if the cardinality of $\mathbb{K}$ is too large, meaning that only a portion of wave vectors can be targeted. Indeed, if the number of independent constraints is larger than the total number of degrees of freedom $dN$, then the system “runs out of degrees of freedom” and the potential energy surface often becomes so complicated that one cannot find an $\Phi=0$ state, even if the target $\mathcal{S}_{0}(\mathbf{k})$ is known to be realizable (Batten et al., 2008). This method also enforces $\mathcal{S}(\mathbf{k})=\mathcal{S}_{0}(\mathbf{k})$ for $\mathbf{k}\in\mathbb{K}$ for a single configuration, while in many cases we only expect the ensemble average $S(\mathbf{k})$ to be equal to $S_{0}(\mathbf{k})$. These drawbacks will be overcome by our new algorithm, as detailed in Sec. II.3.

The discussion above suggests that targeting an ensemble-averaged structure factor, which would enable the toleration of fluctuations in individual configurations, may be a possible way to bypass the limitations of the standard collective-coordinate procedure for the realizability problem. We now show on theoretical grounds how this is indeed the case. Specifically, we demonstrate that targeting ensemble-averaged structure factors results in an enormous increase in the number of degrees of freedom, which in turn enables one to extend the range of constrained wave vectors over an infinite set, in principle. Moreover, we show that the configurations are sampled from a canonical ensemble with a certain pair potential.

Let us begin by imagining imposing constraints such that the average structure factor for a finite number of configurations $N_{c}$ is equal to a target functional form for $\mathbf{k}\in\mathbb{K}$, i.e.,

$$\langle\mathcal{S}(\mathbf{k})\rangle=S_{0}(\mathbf{k})\mbox{, for any }\mathbf{k}\in\mathbb{K},$$

(18)

where $\langle\mathcal{S}(\mathbf{k})\rangle$ is the average structure factor of these $N_{c}$ configurations. We will assume the $N_{c}\to+\infty$ limit in this theoretical subsection, so that

$$\langle\mathcal{S}(\mathbf{k})\rangle=S(\mathbf{k}),$$

(19)

where $S(\mathbf{k})$ is the ensemble-averaged structure factor defined in (2). In this thermodynamic limit ($N_{c}\to+\infty$), there is an infinite number of degrees of freedom, which enables one to extend the range of constrained wave vectors over all space.

A critical question is what is the behavior of each individual configuration when such constraints are imposed? Although it appears that these configurations interact with each other in some complex manner, we now demonstrate that these configurations follow the canonical-ensemble distribution of a pairwise additive potential energy. To begin with, let us consider constraining the target structure factor at a single point, $\mathbf{k}=\mathbf{q}$:

$$S_{0}(\mathbf{q})=\frac{\mathcal{S}_{1}(\mathbf{q})+\mathcal{S}_{2}(\mathbf{q})+\cdots+\mathcal{S}_{N_{c}}(\mathbf{q})}{N_{c}},$$

(20)

where $\mathcal{S}_{i}(\mathbf{q})$ ($i>0$) is the structure factor of the $i$th configuration at wave vector $\mathbf{k}=\mathbf{q}$. In Eq. (20), we are treating $\mathcal{S}_{1}(\mathbf{q})$, $\mathcal{S}_{2}(\mathbf{q})$, $\cdots,\mathcal{S}_{N_{c}}(\mathbf{q})$ as $N_{c}$ random variables and constrain their arithmetic mean to be equal to $S_{0}(\mathbf{q})$. As we will show later, $\mathcal{S}_{i}(\mathbf{q})$ ($i>0$) is exponentially distributed, which implies that $\mathcal{S}_{i}(\mathbf{q})$ is not self-averaging kirkpatric1989 ; parisi2004scale .

The key idea is that because we allow an arbitrary large number of configurations $N_{c}$ and constrain their arithmetic mean structure factor [right-hand side of Eq. (20)], we can focus on one configuration, which we call the “reference system,” with fictitious energy $E_{R}=\mathcal{S}(\mathbf{q})$ and treat the rest of the $N_{c}-1$ configurations as a heat bath with temperature

$$k_{B}T=\frac{S_{0}(\mathbf{q})}{1-S_{0}(\mathbf{q})},$$

(21)

Relation (21) for the temperature is derived immediately below. Under these conditions, the reference system obeys the distribution function of a canonical ensemble in the limit $N_{c}\to\infty$. Figure 1 schematically describes this scenario.

The probability density function (PDF) of the total energy of any system in the canonical ensemble is given by

$$P(E)=\frac{g(E)\exp(-E/k_{B}T)}{Z},$$

(22)

where $Z$ is the partition function, $g(E)$ is energy degeneracy or density of states, and $T$ is the temperature. To determine the temperature of the heat bath explicitly in terms of $S_{0}(\mathbf{q})$, we make the simple observation that in the high-temperature or ideal-gas limit, the PDF (22) is proportional to the density of states

$$g(E)\propto\lim_{T\to\infty}P(E).$$

(23)

Thus, to determine the heat bath’s density of states and corresponding temperature for general correlated systems, we only need to know the distribution of its total energy, $E_{H}$, assuming the heat bath consists of $N_{c}-1$ uncorrelated (infinite-temperature) systems.

Let us first focus on one system comprising the heat bath, and denote its energy by $E_{1}$. For an ideal gas, each particle’s location is random and independent. Thus, for any particle $j$, $\exp(-i\mathbf{q}\cdot\mathbf{r}_{j})$ is a unit vector of a random orientation in the complex plane. Therefore, the collective density variable (17) is the sum of $N$ random unit vectors in the complex plane. From the theory of random walks, we know that for large $N$, ${\tilde{\rho}}(\mathbf{q})$ in the complex plane follows a Gaussian distribution, and the PDF of $|{\tilde{\rho}}(\mathbf{q})|$ is given by

$$P_{uncorr}(|{\tilde{\rho}}(\mathbf{q})|)=\frac{2|{\tilde{\rho}}(\mathbf{q})|}{N}\exp(-|{\tilde{\rho}}(\mathbf{q})|^{2}/N).$$

(24)

Therefore, the PDF of the energy $E_{1}=\mathcal{S}(\mathbf{q})=|{\tilde{\rho}}(\mathbf{q})|^{2}/N$ is given by

$$\begin{split}P_{uncorr}(E_{1})&=P_{uncorr}(|{\tilde{\rho}}(\mathbf{q})|)\left[\frac{d\mathcal{S}(\mathbf{q})}{d|{\tilde{\rho}}(\mathbf{q})|}\right]^{-1}\\ &=\exp(-E_{1}).\end{split}$$

(25)

Thus, the energy (i.e., the single-configuration structure factor) is exponentially distributed. This combined with (23) implies that the density of states of a single configuration is also an exponential function of the energy

$$g(E_{1})\propto\exp(-E_{1}).$$

(26)

For two uncorrelated configurations, the probability distribution of their total energy $E_{12}=E_{1}+E_{2}$ is

	$\displaystyle P_{uncorr}(E_{12})$	$\displaystyle=\int_{0}^{E_{12}}P_{uncorr}(E_{1})P_{uncorr}(E_{2})dE_{1}$		(27)
		$\displaystyle=\int_{0}^{E_{12}}\exp(-E_{1})\exp\left[-(E_{12}-E_{1})\right]dE_{1}$		(28)
		$\displaystyle=E_{12}\exp(-E_{12}).$		(29)

The distribution for the total energy of three configurations is then

$$\begin{split}&P_{uncorr}(E_{123}=E_{12}+E_{3})\\ &=\int_{0}^{E_{123}}P_{uncorr}(E_{3})P_{uncorr}(E_{12})dE_{3}\\ &=\frac{E_{123}^{2}}{2}\exp(-E_{123}).\end{split}$$

(30)

Similarly, one can find that the distribution of the total energy of $N_{c}-1$ configuration is

$$P_{uncorr}(E=E_{1}+E_{2}+\cdots+E_{N_{c}-1})=\frac{E^{N_{c}-1}}{(N_{c}-1)!}\exp(-E).$$

(31)

As we detailed before, the density of states of the heat bath, made from $N_{c}-1$ systems, is proportional to the probability distribution function of the total energy of $N_{c}-1$ uncorrelated systems:

$$g(E)\propto P_{uncorr}(E)=\frac{E^{N_{c}-1}}{(N_{c}-1)!}\exp(-E).$$

(32)

The temperature of the heat bath is therefore

$$k_{B}T=\left[\frac{\partial\ln g(E)}{\partial E}\right]^{-1}=\left[\frac{N_{c}-1}{E}-1\right]^{-1},$$

(33)

where $E=E_{H}$ is the energy of the heat bath. On average, each configuration has an energy of $S_{0}(\mathbf{q})$. Therefore, in the $N_{c}\to\infty$ limit, $E_{H}=(N_{c}-1)S_{0}(\mathbf{q})$ and the heat-bath temperature is explicitly given by

$$k_{B}T=\frac{S_{0}(\mathbf{q})}{1-S_{0}(\mathbf{q})},$$

(34)

which is what we set out to prove.

In the $N_{c}\to\infty$ limit, the heat bath is infinitely large, and we can determine the probability density function of the energy of the reference system, $E_{R}$, which is not included in the heat bath, using the canonical distribution function, i.e.,

	$\displaystyle P(E_{R})$	$\displaystyle=\frac{g(E_{R})\exp(-E_{R}/k_{B}T)}{Z}$		(35)
		$\displaystyle\propto\exp(-E_{R})\exp(-E_{R}/k_{B}T),$		(36)

After normalization, one finds

$$P(E_{R})=\frac{\exp(-E_{R})\exp(-E_{R}/k_{B}T)}{S_{0}(\mathbf{q})}=\frac{\exp[-E_{R}/S_{0}(\mathbf{q})]}{S_{0}(\mathbf{q})}.$$

(37)

Since we previously defined $E_{R}=\mathcal{S}(\mathbf{q})$,

$$P[\mathcal{S}(\mathbf{q})]=\frac{\exp[-\mathcal{S}(\mathbf{q})/S_{0}(\mathbf{q})]}{S_{0}(\mathbf{q})}.$$

(38)

By symmetry, this distribution is applicable not only to the reference system, but also to the other $N_{c}-1$ systems as well. This means that for any particular configuration, its structure factor at a constrained $\mathbf{k}$ vector is exponentially distributed. We will numerically verify this conclusion in the Appendix. As in the ideal-gas case, the exponential distribution of $S(\mathbf{q})$ implies that ${\tilde{\rho}}(\mathbf{q})$ is Gaussian distributed.

As we previously showed, if one constrains the ensemble-averaged structure factor $S(\mathbf{q})$ to be equal to $S_{0}(\mathbf{q})$, the resulting configurations follow canonical-ensemble distribution of a system with energy defined as $E=\mathcal{S}(\mathbf{q})$ at temperature $k_{B}T=S_{0}(\mathbf{q})/(1-S_{0}(\mathbf{q}))$. Equivalently, one could also define a rescaled energy as

$$E={\tilde{v}}(\mathbf{q})\mathcal{S}(\mathbf{q}),$$

(39)

where

$${\tilde{v}}(\mathbf{q})=1/S_{0}(\mathbf{q})-1$$

(40)

and set $k_{B}T={\tilde{v}}(\mathbf{q})S_{0}(\mathbf{q})/(1-S_{0}(\mathbf{q}))=1$. As detailed in our earlier papers Fan et al. (1991); Uche et al. (2004); Torquato et al. (2015), such a definition of $E$ is equivalent to a pairwise additive potential. For example, in the thermodynamic limit, the total energy of such a system of particles with pairwise interactions is given by

$$E=\frac{\rho}{2}\int_{\mathbb{R}^{d}}v(\mathbf{r})g_{2}(\mathbf{r})d\mathbf{r},$$

(41)

which can be represented in Fourier space using Parseval’s theorem Torquato et al. (2015):

$$E=\frac{\rho}{2}{\tilde{v}}(\mathbf{k}=\mathbf{0})-\frac{1}{2}v(\mathbf{r}=\mathbf{0})+\frac{\rho}{2(2\pi)^{d}}\int_{\mathbb{R}^{d}}{\tilde{v}}(\mathbf{k})S(\mathbf{k})d\mathbf{k}.$$

(42)

The derivation of (39) and (40) applies to the cases where one constrains $S(\mathbf{k})$ at a single $\mathbf{k}$ vector. Can one generalize it to constraining $S(\mathbf{k})$ at multiple $\mathbf{k}$ vectors? If one constrains $S(\mathbf{k})$ at up to $d$ different orthogonal wave vectors (inner product being zero), formulas (39) and (40) would still apply exactly. This is because such constraints affect particle positions in different, independent directions, and can thus be treated separately. If the $d$ wave vectors are linearly independent but not orthogonal, one could still apply a linear transformation to reduce the problem to the orthogonal case.

To treat cases where the number of $\mathbf{k}$ vectors is larger than $d$, we recall that in Gibbs formalism, the inverse temperature is a Lagrange multiplier of energy. For multiple constrained wave vectors, $\mathbf{k}_{1}$, $\mathbf{k}_{2}$, $\cdots,\mathbf{k}_{N_{k}}$, we can use a separate Lagrange multiplier for each constraint. Consider maximizing the Gibbs entropy of the reference configuration

$$\mathscr{S}=-k_{B}\int P(\mathbf{r}^{N})\ln P(\mathbf{r}^{N})d\mathbf{r}^{N},$$

(43)

where $P(\mathbf{r}^{N})$ is the probability density function of the reference configuration, subject to constraints

$$C_{\mathbf{k}}=\int P(\mathbf{r}^{N})\mathcal{S}(\mathbf{k})d\mathbf{r}^{N}-S_{0}(\mathbf{k})=0$$

(44)

for each constrained $\mathbf{k}$, and

$$D=\int P(\mathbf{r}^{N})d\mathbf{r}^{N}-1=0,$$

(45)

If these constraints are satisfiable, we can find $P(\mathbf{r}^{N})$ by constructing the Lagrangian function

$$L=\mathscr{S}-\sum_{\mathbf{k}}\lambda_{\mathbf{k}}C_{\mathbf{k}}-\lambda_{D}D,$$

(46)

where $\lambda_{\mathbf{k}}$ and $\lambda_{D}$ are Lagrange multipliers. Setting $\delta L/\delta P(\mathbf{r}^{N})=0$, we find

$$P(\mathbf{r}^{N})=\frac{1}{Z}\exp\left[-\sum_{\mathbf{k}}\lambda_{\mathbf{k}}\mathcal{S}(\mathbf{k})\right],$$

(47)

where $Z=\int\exp\left[-\sum_{\mathbf{k}}\lambda_{\mathbf{k}}\mathcal{S}(\mathbf{k}_{j})\right]d\mathbf{r}^{N}$ is the partition function of the reference system. We see that if we define the fictitious energy as

$$E=\sum_{\mathbf{k}}\lambda_{\mathbf{k}}\mathcal{S}(\mathbf{k}),$$

(48)

which is still a pair potential Fan et al. (1991); Uche et al. (2004); Torquato et al. (2015), then $P(\mathbf{r}^{N})$ follows the equilibrium distribution at $k_{B}T=1$. However, we could not find an explicit expression for $\lambda_{\mathbf{k}}$. Theoretically, $\lambda_{\mathbf{k}}$ is completely determined by the target structure factor. However, the dependency is nontrivial due to the correlation between $S(\mathbf{k})$ at different $\mathbf{k}$ vectors. We proved that $S(\mathbf{k})$ is exponentially distributed in the single-constraint or independent-constraint cases. While we cannot prove that when multiple non-independent wave vectors are constrained, we do provide numerical evidence for such behavior in Appendix A.

To summarize, we have proved that if one constrains the ensemble-averaged structure factor at one or multiple $\mathbf{k}$ vectors, and if the constraints are satisfiable, then the resulting configurations are drawn from the canonical ensemble with a pairwise additive interaction (39) or (48). Thus, when we constrain $S(\mathbf{k})$ for all $\mathbf{k}$ vectors to be equal to a structure factor realized by some $n$-body interactions, our method finds an effective pair interaction that mimics the configurations produced by such $n$-body interactions. For the single-constraint case, we showed that the structure factor for a single configuration at the constrained $\mathbf{k}$ vector is exponentially distributed. For the multiple-constraint case, we will also show strong numerical evidence in Appendix A that the exponential distribution still holds. If the structure factor at the constrained $\mathbf{k}$ vectors are independent from one another, then the interaction is given by Eqs. (39)-(40), and the temperature is $k_{B}T=1$. However, if the structure factor at these $\mathbf{k}$ vectors are correlated, then (40) is inexact. We numerically test and verify these conclusions in Appendix A for specific examples.

Based on this theoretical formalism, we can now straightforwardly devise a new algorithm to construct a canonical ensemble of a finite but large number of configurations, $N_{c}$, targeting a particular functional form for the structure factor. Specifically, we minimize the squared difference for $N_{c}$ configurations but simultaneously, i.e.,

$$\mbox{minimize }\Phi(\mathbf{r}^{\cal{N}})=\sum_{\mathbf{k}\in\mathbb{K}}[\langle\mathcal{S}(\mathbf{k})\rangle-S_{0}(\mathbf{k})]^{2},$$

(49)

where ${\cal N}=NN_{c}$. Thus, compared to the standard collective-coordinate procedure Uche et al. (2006b); Batten et al. (2008); Zachary and Torquato (2011); Zhang et al. (2016), which has available $dN$ number of degrees of freedom, the canonical-ensemble-average generalization, enables us to substantially increase the number of degrees of freedom to $dNN_{c}$. Thus, in practice, the range of wave vectors over which we can constrain the structure factor to have a prescribed functional form can be made to be larger and larger by increasing the number of configurations.

Our algorithm involves minimizing a target function that is a sum over all $\mathbf{k}$ vectors within a wavenumber $K$ from the origin. The number of such $\mathbf{k}$ vectors scales as $K^{d}V$, where $V=N/\rho$ is the volume of the system. For each such $\mathbf{k}$ vector, we need to calculate $N_{c}$ structure factor values, each involves a summation over $N$ particles. Thus, the computational cost scales as $K^{d}VNN_{c}=K^{d}\rho^{-1}N^{2}N_{c}$. As $N$ grows, the computational cost grows quadratically and can become very large. However, the calculations for different $\mathbf{k}$ vectors can be carried out in parallel, and we can thus employ multiple GPUs to overcome the high computational cost. GPUs generally perform single-precision calculations faster than double-precision calculations, but we discovered that double precision is necessary for large system sizes ($N>1000$).

For cases in which $N>2000$, we found that the number of iterations needed to minimize the target function becomes computationally costly. By inspecting the intermediate configurations during minimization, we discovered that $S(\mathbf{k})$ near the origin ($\mathbf{k}=\mathbf{0}$) converges to $S_{0}(\mathbf{k})$ at much slower rate than that of $S(\mathbf{k})$ at other wave vectors. It is reasonable to assume that this slow-up is caused by the fact that changing $S(\mathbf{k})$ at $\mathbf{k}\approx\mathbf{0}$ requires long-range particle motions. To improve the convergence speed for small $k$ when $N>2000$, we introduce a weight of $w(k)=1/k$ in the previous objective function and then carry out the following minimization:

$$\mbox{minimize }\Phi(\mathbf{r}^{\cal{N}})=\sum_{\mathbf{k}\in\mathbb{K}}w(k)[\langle\mathcal{S}(\mathbf{k})\rangle-S_{0}(\mathbf{k})]^{2},$$

(50)

As a proof of concept of these modifications, we were able to generate $N_{c}=100$ configurations, each consisting of $N=20000$ particles, targeting the 1D fermionic target structure factor (58), after five days of computation using four NVIDIA Tesla P100 GPUs, as reported in Sec. III.1. Minimizing the objective function usually requires $\sim 10^{4}$ iterations in 1D but only $10^{2}-10^{3}$ iterations in 2D and 3D. Thus, we can generate higher-dimensional configurations with $N=20000$ particles much faster (about 30 times faster with our hardware).

In subsequent sections, we present results using $N_{c}=100$, which is large but still computationally manageable. Unless otherwise stated, each configuration consists of $N=400$ particles in a linear (1D), square(2D), or cubic (3D) simulation boxes with periodic boundary conditions. As we will show in Fig. 3, $N=400$ is large enough to produce pair statistics indistinguishable from $N=20000$ ones. The pair statistics are averaged over $5000$ configurations to reduce statistical fluctuations. Since the weight $w(k)=1/k$ is necessary only for sufficiently large configurations ($N>2000$), we omit it for simplicity. The set $\mathbb{K}$ contain half of all $\mathbf{k}$ vectors such that $0<|\mathbf{k}|<K$, where $K$ is a constant cutoff. We can omit one half of the $\mathbf{k}$ vectors within the range due to the inversion symmetry of the structure factor: $S(-\mathbf{k})=S(\mathbf{k})$. Unless otherwise stated, we use $K=30$ in 1D and $K=15$ in 2D and 3D. We use the low-storage BFGS algorithm Nocedal (1980); Liu and Nocedal (1989); Johnson to minimize $\Phi$, starting from random initial configurations. After the minimization, $\Phi$ is on the order of $10^{-4}-10^{-6}$. Considering that $\Phi$ is a sum over contributions from $N_{k}=10^{3}-10^{4}$ wave vectors, the difference between $S(\mathbf{k})$ and $S_{0}(\mathbf{k})$ at a particular wave vector is about $\sqrt{\Phi/N_{k}}=10^{-3.5}-10^{-5}$. The efficiency and accuracy of this algorithm is verified by applying to a variety of target structure factors described in Secs. III-VII. We present justifications of this algorithm and parameter choices in Appendix A.

Circular ensembles in the theory of random matrices Mehta (2004) are measures on spaces of unitary matrices introduced by Dyson as modifications of Gaussian matrix ensembles. These different ensembles are equivalent to one another when the matrix size tends to infinity. Dyson showed that the distribution of eigenvalues can be mapped to systems of particles on a circle interacting with a two-dimensional Coulomb potential (logarithmic function) at positive temperatures. These systems in turn can be mapped to logarithmically interacting particles in $\mathbb{R}$ with an appropriately confining potential.

The structure factors that correspond to those of the circular orthogonal ensemble (COE), circular unitary ensemble (CUE), and circular symplectic ensemble (CSE), respectively Mehta (2004); Torquato et al. (2008); Scardicchio et al. (2009); Forrester (2016) at unit density ($\rho=1$) in the thermodynamic limit are:

$\displaystyle S(k)=\left\{\begin{array}[]{lr}{\displaystyle\frac{k}{\pi}-\frac{k}{2\pi}\ln(k/\pi+1)},\quad 0\leq k\leq 2\pi\\ \\ \displaystyle{2-\frac{k}{2\pi}\ln\left(\frac{k/\pi+1}{k/\pi-1}\right)},\quad k>2\pi,\end{array}\right.$

(54)

$\displaystyle S(k)=\left\{\begin{array}[]{lr}{\displaystyle\frac{k}{2\pi}},\quad 0\leq k\leq 2\pi\\ \\ \displaystyle{1},\quad k>2\pi,\end{array}\right.$

(58)

and

$\displaystyle S(k)=\left\{\begin{array}[]{lr}{\displaystyle\frac{k}{4\pi}-\frac{k}{8\pi}\ln\Big{|}1-k/(2\pi)\Big{|}},\quad 0\leq k\leq 4\pi\\ \\ \displaystyle{1},\quad k>4\pi.\end{array}\right.$

(62)

These ensembles correspond to the following values of the inverse temperature $\beta=(k_{B}T)^{-1}$: $\beta=1$ (COE), $\beta=2$ (CUE), and $\beta=4$ (CSE). In all cases, we see that the structure factor $S(k)$ tends to zero linearly in $k$ in the limit $k\to 0$ and hence, according to (10) and (14), are hyperuniform of class II Torquato (2018). The case $\beta=2$ has been generalized to higher dimensions (detailed below) and found to possess identical distribution with a spin-polarized fermionic system Torquato et al. (2008).

The corresponding pair correlations of the COE, CUE and CSE are respectively

$$\begin{split}g_{2}(r)=&1-\frac{\sin^{2}(\pi r)}{(\pi r)^{2}}\\ &+\frac{\Big{(}\pi r\cos(\pi r)-\sin(\pi r)\Big{)}\Big{(}2\,\mbox{Si}(\pi r)-\pi\Big{)}}{2(\pi r)^{2}},\ \end{split}$$

(63)

$\displaystyle g_{2}(r)=1-\frac{\sin^{2}(\pi r)}{(\pi r)^{2}},$

(64)

and

$$\begin{split}g_{2}(r)=&1-\frac{\sin^{2}(2\pi r)}{(2\pi r)^{2}}\\ &+\frac{\Big{(}2\pi r\cos(2\pi r)-\sin(2\pi r)\Big{)}\,\mbox{Si}(2\pi r)}{4(\pi r)^{2}},\end{split}$$

(65)

We now apply our algorithm by targeting these three structure factors for systems with $N=400$. The target analytical forms for the structure factor for $\beta=1$, $\beta=2$, and $\beta=4$ and the corresponding simulation data are plotted in Figs.  2,  3,  and  4, respectively. We include in these figures the corresponding pair correlation functions, both the analytical forms and the simulation data, as obtained by sampling the generated configurations. From all of these figures, we see that our targeting algorithm enables us to realize these ensembles with high accuracy, validating its utility and applicability. For the $\beta=2$ case, we also carried out the simulation results for a much larger system with $N=20,000$. We see from Fig. 3 that the results for $N=400$ are indistinguishable from those for $N=20000$.

By a theorem of Henderson Henderson (1974), a pair potential that gives rise to a given pair correlation function is unique up to a constant shift, although this cannot apply at $T=0$ foo ; Cohn and Kumar (2007); Cohn et al. (2019). The fact that our methodology yields ensembles of configurations with targeted pair statistics (whenever realizable) that are determined by effective pair potentials means that those interactions in the case of Dyson’s one-dimensional COE, CUE and CSE must exactly be given by the two-dimensional Coulombic potential.

Costin and Lebowitz Costin and Lebowitz (2004) showed that there exists a one-dimensional determinantal point process at unit density in which the total correlation function is the following simple exponential function:

$$h(r)=-\exp(-\lambda r),$$

(66)

where $\lambda\geq 2$. This corresponds to a structure factor with a Lorentzian form:

$$S(k)=\frac{\lambda(\lambda-2)+k^{2}}{k^{2}+\lambda^{2}}.$$

(67)

This result implies that the system is hyperuniform at the “borderline” case of $\lambda=2$, since the structure factor $S(k)=k^{2}/(4+k^{2})$ tends to zero quadratically in $k$ in the limit $k\to 0$ and hence, according to (10) and (14), are hyperuniform of class I Torquato (2018). We have targeted the $\lambda=2$ target, and successfully realize it. The results for the pair correlation function and the structure factor are presented in Fig. 5.

The 1D CUE point process with a structure factor given by Eq. (58) has been generalized to so-called “Fermi-sphere” point processes in $d$-dimensional Euclidean space $\mathbb{R}^{d}$ Torquato et al. (2008). Specifically, such disordered hyperuniform point processes correspond to the spatial distribution of spin-polarized free fermions in $\mathbb{R}^{d}$, which are special cases of determinantal processes. In particular, the structure factor in $\mathbb{R}$ at unit density is given by

$$S(k)=1-\alpha(k,\kappa),$$

(68)

where $\alpha(k,\kappa)$ is the volume common to two spherical windows of radius $\kappa$ whose centers are separated by a distance $k$ divided by $v_{1}(\kappa)$, the volume of a spherical window of radius $\kappa=2\sqrt{\pi}[\Gamma(1+d/2)]^{1/d}$, which is known analytically in any dimension Torquato and Stillinger (2006a). This result implies that the structure factor $S(k)$ tends to zero linearly in $k$ in the limit $k\to 0$ and hence are hyperuniform of class II Torquato (2018). The corresponding pair correlation function of such a point process is given by

$\displaystyle g_{2}(r)=1-2^{d}\Gamma(1+d/2)^{2}\frac{J^{2}_{d/2}(\kappa r)}{(\kappa r)^{d}},$

(69)

where $J_{\nu}(x)$ is the Bessel function of the first kind of order $\nu$.

We have applied our algorithm to target the structure factor (68) in two and three dimensions using $K=15$. The results are presented in Figs. 6 and 7 along with the corresponding pair correlation functions sampled from the generated configurations as well as the analytical forms obtained from (69). Consistent with the known realizability of these targets, we see excellent agreement between the targeted structure factors and those obtained from our ensemble-average formulation. A two-dimensional configuration is shown in Fig. 8.

Unlike the one-dimensional COE, CUE, and CSE determinantal point configurations, the interaction potential for general determinantal point processes must contain at least up to three-body potentials; see the Appendix of Ref. Torquato et al., 2008. Thus, for the 2D and 3D fermi-sphere targets as well as the 1D Lorentzian target (Sec. III.2), we show for the first time that there exists effective pair interactions that mimic the higher-order $n$-body interactions corresponding to these determinantal point processes.

An example of a 2D determinantal point process that exhibits hyperuniform behavior is generated by the Ginibre ensemble Jancovici (1981), which is a special case of the two-dimensional one-component plasma Jancovici (1981). A one-component plasma (OCP) is an equilibrium system of identical point particles of charge $e$ interacting via the log Coulomb potential and immersed in a rigid, uniform background of opposite charge to ensure overall charge neutrality. For $\beta=2$, the total correlation function for the OCP (Ginibre ensemble) in the thermodynamic limit was found exactly by Jancovici Jancovici (1981):

$\displaystyle h(r)=-\exp\left(-\rho\pi r^{2}\right).$

(70)

The corresponding structure factor is given by

$$S(k)=1-\exp\left(-\frac{k^{2}}{4\pi\rho}\right).$$

(71)

This result implies that the structure factor $S(k)$ tends to zero quadratically in $k$ in the limit $k\to 0$ and hence are hyperuniform of class II Torquato (2018).

Using our method, we targeted the OCP structure factor (71) using $K=15$. As shown in Fig. 9, it is seen that the algorithm is able to realize this target with very high accuracy. The corresponding pair correlation function obtained by sampling the resulting configurations agrees very well with the exact $g_{2}(r)$ obtained from (70), as shown in Fig. 9. One configuration is shown in Fig. 10. It it noteworthy that we had previously employed a completely different algorithm to generate these configurations as well as other determinantal point processes Scardicchio et al. (2009), but the maximum attainable system sizes were substantially much smaller ($N\approx 100$) in that study.

To begin, we ask whether a total correlation function with the following Gaussian form is realizable as a hyperuniform system across dimensions:

$$h(r)=-\exp\left(-(r/a)^{2}\right),$$

(74)

where $a$ is a positive constant. The corresponding structure factor is given by

$$S(k)=1-\rho a^{d}\pi^{d/2}\exp\left(-\frac{k^{2}}{4\pi\rho}\right),$$

(75)

which implies that $S(k)$ tends to zero quadratically as $k\to 0$ for all $d$. The hyperuniformity condition requires the unique density to be given by $\rho=(a\sqrt{\pi})^{-d}$. Note that the case $d=2$ is exactly the same as the two-dimensional OCP system in which $h(r)$ and $S(k)$ are given by (70) and (71), respectively.

We target such structure factors in one and three dimensions and find that they are realizable as hyperuniform systems at unit density. For $d=1$, we find excellent agreement between the simulated and target structure factors are obtained, as shown in Fig. 12. This strongly suggests that such systems are realizable in one dimension, the most difficult dimensionality case. Indeed, we also find the same excellent agreement between the simulated and target structure factors in three dimensions, as illustrated in Fig. 13. We conclude that such targets are realizable as disordered hyperuniform systems of class I Torquato (2018) in any space dimension whenever $\rho=(a\sqrt{\pi})^{-d}$. A 3D configuration is shown in Fig. 14.

Consider the following $d$-dimensional generalization of the total correlation function of the OCP:

$$h(r)=-\exp[-\rho v_{1}(r)],$$

(76)

where $v_{1}(r)$ is the volume of a sphere of radius $r$ [cf. (4)]. It is noteworthy that such a total correlation function automatically satisfies the hyperuniformity requirement for any positive density and any $d$, since ${\tilde{h}}(k=0)=\int_{\mathbb{R^{d}}}h(r)d{\bf r}=-1/\rho$ [cf. (2) and (5)]. Note that when $d=1$, this is identical to the realizable total correlation function (66) with $\lambda=2\rho$. Moreover, when $d=2$, this is identical to the realizable OCP function (70).

It is not known whether configurations corresponding to (76) for $d\geq 3$ are realizable. We target the structure factor in this case in three dimensions:

$$\begin{split}S(k)=&1-\frac{1}{1080\rho^{4/3}\pi^{2/3}\Gamma(2/3)}\bigg{\{}\\ &2^{1/3}3^{5/6}\pi^{1/3}{}_{0}F_{3}\left[-;\frac{4}{3},\frac{3}{2},\frac{11}{6};a(k)\right]k^{4}\\ &-30(6\rho)^{2/3}\Gamma\left(\frac{2}{3}\right)^{2}{}_{0}F_{3}\left[-;\frac{2}{3},\frac{7}{6},\frac{3}{2};a(k)\right]k^{2}\\ &+1080\Gamma\left(\frac{2}{3}\right)\pi^{2/3}\rho^{4/3}{}_{1}F_{4}\left[1;\frac{1}{3},\frac{2}{3},\frac{5}{6},\frac{7}{6};a(k)\right]\bigg{\}},\end{split}$$

(77)

where $\Gamma(x)$ is the gamma function, ${}_{p}F_{q}(a_{1},\cdots,a_{p};b_{1},\cdots,b_{q};z)$ is the generalized hypergeometric function, and

$$a(k)=\frac{k^{6}}{20736\pi^{2}\rho^{2}}.$$

(78)

For small wavenumbers, this 3D structure factor as well as those for any other values of $d$ goes to zero quadratically in $k$ as $k$ tends to zero; specifically,

$$S(k)\sim k^{2}\qquad(k\to 0).$$

(79)

This means that $S(k)$ is analytic at the origin, which in turn implies that $h(r)$ decays to zero exponentially fast or faster Torquato (2018). We find that this 3D structure factor is indeed realizable for $\rho=1$. The results depicted in Fig. 15 show excellent agreement between the simulated and target structure factors. One such configuration is shown in Fig. 16. Since (76) is realizable for $d=3$, it should be realizable in higher dimensions and hence such systems in $\mathbb{R}^{d}$ for any $d$ are hyperuniform of class I Torquato (2018).